Chapter 4 

Introduction to TDDFT 

Eberhard K. U. Gross and Neepa T. Maitra 



4.1 Introduction 

Correlated electron motion plays a significant role in the spectra described in the 
previous chapters. Further, placing an atom, molecule or solid in a strong laser field 
reveals fascinating non-perturbative phenomena, such as non- sequential multiple- 
ionization (see Chap. 18), whose origins lie in the subtle ways electrons interact with 
each other. The direct approach to treat these problems is to solve the (non-relativistic) 
time-dependent Schrodinger equation for the many-electron wavefunction V(t): 

H(t)^(t)=i—^- 9 H(t) = T + V ee + y ext (0 (4.1) 
ot 

for a given initial wavefunction ^ (0) . Here, the kinetic energy and electron-electron 
repulsion, are, respectively: 

1 N 1*1 
t = --^l and £ , (4 ' 2) 

i=l i*/ 1 1 Jl 

and the "external potential" represents the potential the electrons experience due to 
the nuclear attraction and due to any field applied to the system (e.g. laser): 
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N 



(4.3) 



i=\ 



For example, u e xt(f * , t) can represent the Coulomb interaction of the electrons with 
a set of nuclei, possibly moving along some classical path, 



where Z v and R v denote the charge and position of the nucleus v, and N n stands 
for the total number of nuclei in the system. This may be useful to study, e.g., 
scattering experiments, chemical reactions, etc. Another example is the interaction 
with external fields, e.g. for a system illuminated by a laser beam we can write, in 
the dipole approximation, 



where a, co and E are the polarization, the frequency and the amplitude of the laser, 
respectively. The function f(t) is an envelope that describes the temporal shape of 
the laser pulse. We use atomic units (e 2 = H = m = 1) throughout this chapter; all 
distances are in Bohr, energies in Hartrees (1 H = 27.21 eV = 627.5 kcal/mol), and 
times in units of 2.419 x 10 -17 s. 

Solving Eq. 4.1 is an exceedingly difficult task. Even putting aside 
time-dependence, the problem of finding the ground-state scales exponentially with 
the number of electrons. Moreover, ^ contains far more information than one could 
possibly need or even want. For example, consider storing the ground state of the 
oxygen atom, and for simplicity, disregard spin. Then & depends on 24 coordinates, 
three for each of the eight electrons. Allowing ourselves a modest ten grid-points for 
each coordinate, means that we need 10 24 numbers to represent the wavefunction. 
Assuming each number requires one byte to store, and that the capacity of a DVD is 
10 10 bytes, we see that 10 14 DVD's are required to store just the ground-state wave- 
function of the oxygen atom, even on a coarse grid. Physically, we are instead inter- 
ested in integrated quantities, such as one- or two-body probability-densities, which, 
traditionally can be extracted from this foreboding However, a method that could 
yield such quantities directly, by-passing the need to calculate would be highly 
attractive. This is the idea of density-functional theories. In fact, in 1964, Hohen- 
berg and Kohn (1964), proved that all observable properties of a static many-electron 
system can be extracted exactly, in principle, from the one-body ground-state density 
alone. Twenty years later, Runge and Gross extended this to time-dependent systems, 
showing that all observable properties of a many-electron system, beginning in a 
given initial state & (0) , may be extracted from the one-body time-dependent density 
alone (Runge and Gross 1984). What has made (TD)DFT so incredibly successful is 
the Kohn-Sham (KS) system: the density of the interacting many-electron system is 





(4.5) 
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obtained as the density of an auxiliary system of non-interacting fermions, living in 
a one-body potential. The exponential scaling with system- size that the solution to 
Eq. 4.1 requires is replaced in TDDFT by the much gentler N 3 or N 2 scaling 
(depending on the implementation) (Marques 2006), opening the door to the quantum 
mechanical study of much larger systems, from nanoscale devices to biomolecules. 
(See Chaps. 19-21 for details on the numerical issues). Although the ground-state 
and time-dependent theories have a similar flavor, and modus operandi, their proofs 
and functionals are quite distinct. 

Before we delve into the details of the fundamental theorems of TDDFT in the 
next section, we make some historical notes. As early as 1933, Bloch proposed a time- 
dependent Thomas-Fermi model (Bloch 1933). Ando (1977a, b), Peuckert (1978) 
and Zangwill and Soven (1980a, 1980b) ran the first time-dependent KS calculations, 
assuming a TDKS theorem exists. They treated the linear density response to a 
time-dependent external potential as the response of non-interacting electrons to an 
effective time-dependent potential. Ando calculated resonance energy and absorption 
lineshapes for intersubband transitions on the surface of silicon, while Peuckert and 
Zangwill and Soven studied rare-gas atoms. In analogy to ground- state KS theory, this 
effective potential was assumed to contain an exchange-correlation part, v xc (r, t), 
in addition to the time-dependent external and Hartree terms. Peuckert suggested an 
iterative scheme for the calculation of v xc , while Ando, Zangwill and Soven adopted 
the functional form of the static exchange-correlation potential in LDA. Significant 
steps towards a rigorous foundation of TDDFT were taken by Deb and Ghosh (1982), 
Ghosh and Deb (1982, 1983a, 1983b) for time-periodic potentials and by Bartolotti 
(1981, 1982, 1984, 1987) for adiabatic processes. These authors formulated and 
explored Hohenberg-Kohn and KS type theorems for the time-dependent density 
for these cases. Modern TDDFT is based on the general formulation of Runge and 
Gross (1984). 

TDDFT is being used today in an ever-increasing range of applications to widely- 
varying systems in chemistry, biology, solid-state physics, and materials science. 
We end the introduction by classifying these into four areas. First, the vast majority of 
TDDFT calculations today lie in spectroscopy (Chaps. 1-3 earlier), yielding response 
and excitations of atoms, molecules and solids. The laser field applied to the system, 
initially in its ground-state, is weak and perturbation theory applies. We need to know 
only the exchange-correlation potential in the vicinity of the ground-state, and often 
(but not always), formulations directly in frequency-domain are used (Chap. 7 and 
later in this chapter). Overall, results for excitation energies tend to be fairly good 
(few tenths of an eV error, typically) but depend significantly on the system and type 
of excitation considered, e.g. the errors for long-range charge-transfer excitations 
can be ten times as large. For solids, to obtain accurate optical absorption spectra 
of insulators needs functionals more sophisticated than the simplest ones (ALDA), 
however ALDA does very well for electron-energy-loss spectra. We shall return to 
the general performance of TDDFT for spectra in Sect. 4.8. 

The second class of applications is real-time dynamics in non-perturbative fields. 
The applied electric field is comparable to, or greater than, the static electric field due 
to the nuclei. Fascinating and subtle electron interaction effects can make the "single 
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active electron" picture often used for these problems break down. For dynamics 
in strong fields, it is pushing today's computational limits for correlated wavefunc- 
tion methods to go beyond one or two electrons in three-dimensions, so TDDFT is 
particularly promising in this regime. However, the demands on the functionals 
for accurate results can be challenging. Chapter 18 discusses many interesting 
phenomena, that reveal fundamental properties of atoms and molecules. Also under 
this umbrella are coupled electron-ion dynamics. For example, applying a laser 
pulse to a molecule or cluster, drives both the electronic and nuclear system out 
of equilibrium; generally their coupled motion is highly complicated, and various 
approximation schemes to account for electronic-nuclear "back-reaction" have been 
devised. Often in photochemical applications, the dynamics are treated beginning 
on an excited potential energy surface; that is, the dynamics leading to the initial 
electronic excitation is not explicitly treated but instead defines the initial state for 
the subsequent field- free dynamics of the full correlated electron and ion system. We 
refer the reader to Chaps. 14-16 for both formal and practical discussions on how to 
treat this challenging problem. 

The third class of applications returns to the ground-state: based on the fluctuation- 
dissipation theorem, one can obtain an expression for the ground-state exchange- 
correlation energy from a TDDFT response function, as is discussed in Chap. 22. Such 
calculations are significantly more computationally demanding than usual ground- 
state calculations but provide a natural methodology for some of the most difficult 
challenges for ground-state approximations, in particular van der Waal's forces. We 
refer the reader to Chaps. 22 and 23 for a detailed discussion. 

The fourth class of applications is related to viscous forces arising from electron- 
electron interactions in very large finite systems, or extended systems such as solids. 
Consider an initial non-equilibrium state in such a system, created, for example, by 
a laser pulse that is then turned off. For a large enough system, electron interac- 
tion subsequently relaxes the system to the ground-state, or to thermal equilibrium. 
Relaxation induced by electron-interaction can be in principle exactly captured in 
TDDFT, but for a theoretically consistent formulation, one should go beyond the 
most commonly used approximations. A closely related approach is time-dependent 
current-density functional theory (TDCDFT) (Sect. 4.4.4 and Chap. 24), where an 
xc vector-potential provides the viscous force (D' Agosta and Vignale 2006; Ullrich 
2006b). Dissipation phenomena studied so far, using either TDDFT or TDCDFT, 
include energy loss in atomic collisions with metal clusters (Baer and Siam 2004), 
the stopping power of ions in electron liquids (Nazarov et al. 2005; Hatcher 2008; 
Nazarov 2007), spin-Coulomb drag (d'Amico and Ullaich 2006, 2010), and a hot 
electron probing a molecular resonance at a surface (Gavnholt et al. 2009). Transport 
through molecular devices (Chap. 17) is an important subgroup of these applications: 
how a system evolves to a steady-state after a bias is applied (Stefanucci Almbhdh 
2004a; Kurth et al. 2005, 2010; Khosravi et al. 2009; Stefanucci 2007; Koentopp 
et al. 2008; Zheng 2010). In this problem, to account for coupling of electrons in 
the molecular wire to a "bath" of electrons in the leads, or to account for coupling 
to external phonon modes, one is led to the "open systems" analyses, reviewed in 
Chaps. 10 and 11 (Gebauer and Car 2004a; Burke 2005; Yuen-Zhou et al. 2010; 
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Di Ventra D' Agosta 2007; Chen 2007; Appel and di ventra 2009). In Ullrich (2002b) 
the linear response of weakly disordered systems was formulated, embracing both 
extrinsic damping (interface roughness and charged impurities) as well as intrinsic 
dissipation from electron interaction. 

The present chapter is organized as follows. Section 4.2 presents the proof of the 
fundamental theorem in TDDFT. Section 4.3 then presents the time-dependent KS 
equations, the "performers" of TDDFT. In Sect. 4.4 we discuss several details of 
the theory, which are somewhat technical but important if one scratches below the 
surface. Section 4.5 derives the linear response formulation, and the matrix equations 
which run the show for most of the applications today. Section 4.6 briefly presents the 
equations for higher-order response, while Sect. 4.7 describes some of the approxi- 
mations in use currently for the xc functional. Finally, Sect. 4.8 gives an overview of 
the performance and challenges for the approximations today. 



4.2 One-to-One Density-Potential Mapping 

The central theorem of TDDFT (the Runge-Gross theorem) proves that there is a one- 
to-one correspondence between the external (time-dependent) potential, v QXt (r, t), 
and the electronic one-body density, n(r, t), for many-body systems evolving from 
a fixed initial state (Runge and Gross 1984). The density n(r,t) is the probability 
(normalized to the particle number AO of finding any one electron, of any spin cr, at 
position r : 



The implications of this theorem are enormous: if we know only the time-dependent 
density of a system, evolving from a given initial state, then this identifies the external 
potential that produced this density. The external potential completely identifies the 
Hamiltonian (the other terms given by Eq. 4.2 are determined from the fact that 
we are dealing with electrons, with N being the integral of the density of Eq. 4.6 
over r.) The time-dependent Schrodinger equation can then be solved, in principle, 
and all properties of the system obtained. That is, for this given initial- state, the 
electronic density, a function of just three spatial variables and time, determines all 
other properties of the interacting many-electron system. 

This remarkable statement is the analogue of the Hohenberg-Kohn theorem for 
ground- state DFT, where the situation is somewhat simpler: the density -potential 
map there holds only for the ground-state, so there is no time-dependence and no 
dependence on the initial state. The Hohenberg-Kohn proof is based on the Rayleigh- 
Ritz minimum principle for the energy. A straightforward extension to the time- 
dependent domain is not possible since a minimum principle is not available in this 
case. 




(4.6) 
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Instead, the proof for a 1-1 mapping between time-dependent potentials and 
time-dependent densities is based on considering the quantum-mechanical equation 
of motion for the current-density, for a Hamiltonian of the form of Eqs. 4.1-4.3. 
The proof requires the potentials v QXt (r, t) to be time-analytic around the initial 
time, i.e. that they equal their Taylor-series expansions in t around t = 0, for a finite 
time interval: 

OO j 

V«t(r, = X 7^exU(r)^. (4.7) 
k=0 

The aim is to show that two densities n(r,t) and n f (r, t) evolving from a common 
initial state under the influence of the potentials v ext (r,t) and v f ext (r,t) are always 
different provided that the potentials differ by more than a purely time-dependent 
function: 

v ext (r,t)^v f ext (r,t) + c(t). (4.8) 

The above condition is a physical one, representing simply a gauge-freedom. A purely 
time-dependent constant in the potential cannot alter the physics: if two potentials 
differ only by a purely time-dependent function, their resulting wavefunctions differ 
only by a purely time-dependent phase factor. Their resulting densities are iden- 
tical. All variables that correspond to expectation values of Hermitian operators are 
unaffected by such a purely time-dependent phase. There is an analogous condi- 
tion in the ground-state proof of Hohenberg and Kohn. Equation 4.8 is equivalent 
to the statement that for the expansion coefficients u e xt,fc(f ) an d v' QXt k (r) [where, 

as in Eq. 4.7, v f exi (r, = Xi£o H^ext fc( r )** 1 there exists a smallest integer k > 
such that 



VextA r ) - Kxt,k( r ) = [^xt(f, t) - v' ext (r, t)] 



^ const. (4.9) 

=o 



The initial state need not be the ground-state or any stationary state of the initial 
potential, which means that "sudden switching" is covered by the RG theorem. But 
potentials that turn on adiabatically from t = — oo, are not, since they do not satisfy 
Eq. 4.7 (see also Sect. 4.5.3). 

The first step of the proof demonstrates that the current-densities 

j(r,t) = (V(t)\j(r)\V(t)) (4.10) 

and 

j'(r,t) = (V'(t)\j(r)\V'(t)) (4.11) 
are different for different potentials v txt and v' GXt . Here, 
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1 N 

j(r) = -X [V/5(r - n) + 5(r - r,-)V/] (4.12) 
i=i 

is the usual paramagnetic current-density operator. In the second step, use of the 
continuity equation shows that the densities n and n' are different. We now proceed 
with the details. 

Step 1 We apply the equation of motion for the expectation value of a general operator 
G(0, 



(tf(f)l Q(t)\V(t)) = (V(t)\ ^ - i[g(f), » (01^ 1^(0), (4.13) 



to the current densities: 

9 9 

-j(r, f) = -(^(f)|}(r)|^(f)> = -imt)\U(r), H (t)]\<P (t)) (4.14a) 
df ot 

9 9 

-/(r, f) = -r-{V'(f)\j(r)\V'(t)) = -i{V'(t)\[j(r), H'(f)W(t)), (4.14b) 
ot ot 

and take their difference evaluated at the initial time. Since V and ^' evolve from 
the same initial state 

\p(t = 0) = V\t = 0) = Vo, (4.15) 

and the corresponding Hamiltonians differ only in their external potentials, we have 
9 



H V(r,t)-f(r,t)] 



= -i(Vo\U(r),H(0)-H f (0)Wo) 
= - n (r)W [ v&A (r, 0) - < xt (r, 0)] (4.16) 



where no(r) is the initial density. Now, if the condition (4.9) is satisfied for k = the 
right-hand side of (4.16) cannot vanish identically and j and ]' will become different 
infinitesimally later than t = 0. If the smallest integer k for which Eq. 4.9 holds is 
greater than zero, we use Eq. 4.13 (k + 1) times. That is, as for k = above where 
we used Q(t) = j(r) in Eq. 4.13, for k = 1, we take Q(t) = -i[j(r), H(t)]; for 
general k, Q(t) = (-i) h [[[j (r), H(t)l H(t)] . . . H(t)] k meaning there are k nested 
commutators to take. After some algebra 1 : 



Jt) Wr,t)-j'(r,m 



= -n (r)Ww k (r)^0 (4.17) 



Note that Eq. 4.17 applies for all integers from to this smallest k for which Eq. 4.9 holds, but 
not for integers larger than this smallest k. 



60 



E. K. U. Gross and N. T. Maitra 



with 



u>k(r) = (^\ [v ext (r, t) - v ext (r, t)] 



(4.18) 



Once again, we conclude that infinitesimally later than the initial time, 

j(r,t)^f(r,t). (4.19) 

This first step thus proves that the current-densities evolving from the same initial 
state in two physically distinct potentials, will differ. That is, it proves a one-to-one 
correspondence between current-densities and potentials, for a given initial- state. 

Step 2 To prove the corresponding statement for the densities we use the continuity 
equation 

dn(r,t) 

—y-L = -V.j(r,t) (4.20) 
ot 

to calculate the (k + 2)nd time-derivative of the density n(r,t) and likewise of the 
density n f (r,t). Taking the difference of the two at the initial time t = and inserting 
Eq. 4.18 yields 



[n(r, t) - n'(r, t)] 



= V ■ [n (r)Vw k (r)l (4.21) 



Now, if there was no divergence-operator on the r.h.s., our task would be complete, 
showing that the densities n(r,t) and n'(r, t) will become different infinitesimally 
later than t = 0. To show that the divergence does not render the r.h.s. zero, thus 
allowing an escape from this conclusion, we consider the integral 



/ 



d 3 rn (r)[Vw k (r)] 2 = - / d 3 rw k (r)V ■ [n (r)Vw k (r)] 



+ <bdS>[no(r)w k (r)Vw k (r)], (4.22) 



where we have used Green's theorem. For physically reasonable potentials (i.e. poten- 
tials arising from normalizable external charge densities), the surface integral on the 
right vanishes (Gross and Kohn 1990) (more details are given in Sect. 4.4.1). Since 
the integrand on the left-hand side is strictly positive or zero, the first term on the 
right must be strictly positive. That is, V • [no(r)Vwk(r)] cannot be zero everywhere. 
This completes the proof of the theorem. 

We have shown that densities evolving from the same initial wavefunction &o 
in different potentials must be different. Schematically, the Runge-Gross theorem 
shows 
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*o : uext «-* n. (4.23) 

The backward arrow, that a given time-dependent density points to a single time- 
dependent potential for a given initial state, has been proven above. The forward arrow 
follows directly from the uniqueness of solutions to the time-dependent Schrodinger 
equation. 

Due to the one-to-one correspondence, for a given initial state, the time-dependent 
density determines the potential up to a purely time-dependent constant. The wave- 
function is therefore determined up to a purely-time-dependent phase, as discussed 
at the beginning of this section, and so can be regarded as a functional of the density 
and initial state: 

qs(t) = c~ ia{t) ^[n, V ](t). (4.24) 

As a consequence, the expectation value of any quantum mechanical Hermitian oper- 
ator Q(t) is a unique functional of the density and initial state (and, not surprisingly, 
the ambiguity in the phase cancels out): 

Q[n, m(t) = (V[n, m(t)\Q(t)\V[n, #b](0>- (4.25) 

We also note that the particular form of the Coulomb interaction did not enter 
into the proof. In fact, the proof applies not just to electrons, but to any system of 
identical particles, interacting with any (but fixed) particle-interaction, and obeying 
either fermionic or bosonic statistics. 

In Sect. 4.4 we shall return to some details and extensions of the proof, but now 
we proceed with how TDDFT operates in practice: the time-dependent Kohn-Sham 
equations. 



4.3 Time-Dependent Kohn-Sham Equations 

Finding functionals directly in terms of the density can be rather difficult. In partic- 
ular, it is not known how to write the kinetic energy as an explicit functional of the 
density. The same problem occurs in ground-state DFT, where the search for accu- 
rate kinetic-energy density-functionals is an active research area. Instead, like in the 
ground- state theory, we turn to a non-interacting system of fermions called the Kohn- 
Sham (KS) system, defined such that it exactly reproduces the density of the true 
interacting system. A large part of the kinetic energy of the true system is obtained 
directly as an orbital-functional, evaluating the usual kinetic energy operator on the 
KS orbitals. (The rest, along with other many-electron effects, is contained in the 
exchange-correlation potential.) All properties of the true system can be extracted 
from the density of the KS system. 

Because the 1-1 correspondence between time-dependent densities and time- 
dependent potentials can be established for any given interaction V ee , in particular 
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also for Vee = 0, it applies to the KS system. Therefore the external potential 
^ksD*; ^o](7% of a non-interacting system reproducing a given density n(r,t), 
starting in the initial state @o, is uniquely determined. The initial KS state is 
almost always chosen to be a single Slater determinant of single-particle spin-orbitals 
(Pi (r , 0) (but need not be) ; the only condition on its choice is that it must be compatible 
with the given density. That is, it must reproduce the initial density and also its first 
time-derivative (from Eq. 4.20, see also Eqs. 4.29-4.30 shortly). However, the 1-1 
correspondence only ensures the uniqueness of uksD*; ^o] but not its existence for 
an arbitrary n(r,t). That is, the proof does not tell us whether a KS system exists 
or not; this is called the non-interacting u-representability problem, similar to the 
ground-state case. We return to this question later in Sect. 4.4.2, but for now we 
assume that uks exists for the time-dependent density of the interacting system of 
interest. Under this assumption, the density of the interacting system can be obtained 
from 

N 

n(r,t) = Y J Wj(r,t)\ 1 (4.26) 

7=1 



with orbitals (Pj(r, t) satisfying the time-dependent KS equation 

V 2 



, vks[«; <t>o\(r,t) 



<Pj(r,t). (4.27) 



dV w(r ' *\ +v xc [n;Vo,$o](r,t), 



Analogously to the ground-state case, vks is decomposed into three terms: 

/ 

(4.28) 

where v ext [n; &o](r,t) is the external time-dependent field. The second term on 
the right-hand side of Eq. 4.28 is the time-dependent Hartree potential, describing 
the interaction of classical electronic charge distributions, while the third term is 
the exchange-correlation (xc) potential which, in practice, has to be approximated. 
Equation 4.28 defines the xc potential: it, added to the classical Hartree potential, 
is the difference between the external potential that generates density n(r, t) in an 
interacting system starting in initial state and the one-body potential that generates 
this same density in a non-interacting system starting in initial state <£>o- 

The functional-dependence of t> e xt displayed in the first term on the r.h.s. of 
Eq. 4.28 is not important in practice, since for real calculations, the external potential 
is given by the physics at hand. Only the xc potential needs to be approximated in 
practice, as a functional of the density, the true initial state and the KS initial state. 
This functional is a very complex one: knowing it implies the solution of all time- 
dependent Coulomb interacting problems. 

As in the static case, the great advantage of the time-dependent KS scheme lies 
in its computational simplicity compared to other methods such as time-dependent 
configuration interaction (Errea et al. 1985; Reading and Ford 1987; Krause 2005; 
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Krause 2007) or multi-configuration time-dependent Hartree-Fock (Zanghellini et 
al. 2003; Kato and Kono 2004; Nest 2005; Meyer 1990; Caillat 2005). A TDKS 
calculation proceeds as follows. An initial set of Af orthonormal KS orbitals is chosen, 
which must reproduce the exact density of the true initial state &o (given by the 
problem) and its first time-derivative: 

rc(r,0) = ^Mr,0)| 2 = Af I d3r 2'" / d 3 r N \V (x,x 2 ...x N )\ 2 

i=l o,o 2 ...o N J J 

(4.29) 

(and we note there exist infinitely many Slater determinants that reproduce a given 
density (Harriman 1981; Zumbach and Maschke 1983)), and 

N 

n(r, 0) = -V • 3m^2>*(r, 0)V^-(r, 0) 

i=\ a 

= -NW -3m ^ / d 3 r 2 --- / d 3 r N ^ (x , x 2 . . . x N )W^ (x , x 2 . . . x N ) 

(4.30) 

using notation x = (r, a). The TDKS equations (Eq. 4.27) then propagate these 
initial orbitals, under the external potential given by the problem at hand, together 
with the Hartree potential and an approximation for the xc potential in Eq. 4.28. In 
Sect. 4.7 we shall discuss the approximations that are usually used here. 

The choice of the KS initial state, and the fact that the KS potential depends on 
this choice is a completely new feature of TDDFT without a ground-state analogue. 
A discussion of the subtleties arising from initial- state dependence can be found 
in Chap. 8. In practice, the theory would be much simpler if we could deal with 
functionals of the density alone. For a large class of systems, namely those where 
both and 0o are non-degenerate ground states, observables are indeed functionals 
of the density alone. This is because any non-degenerate ground state is a unique 
functional of its density no(r) by virtue of the traditional Hohenberg-Kohn theorem 
(Hohenberg and Kohn 1964). In particular, the initial KS orbitals are uniquely deter- 
mined as well in this case. We emphasize this is often the case in practice; in particular 
in the linear response regime, where spectra are calculated (see Sect. 4.5). This is 
where the vast majority of applications of TDDFT lie today. 



4.4 More Details and Extensions 

We now discuss in more detail some important points that arose in the derivation of 
the proof above. Several of these are discussed at further length in the subsequent 
chapters. 
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4.4.1 The Surface Condition 



It is essential for Step 2 of the RG proof that the surface term 



dS • no(r)w k (r)Ww k (r) 



(4.31) 



appearing in Eq. 4.22 vanishes. Let us consider realistic physical potentials of the 
form 



where n QXt (r,t) denotes normalizable charge-densities external to the electronic 
system. A Taylor expansion in time of this expression shows that the coefficients 
v ex t,k, an d therefore wu fall off at least as 1/r asymptotically so that, for physical 
initial densities, the surface integral vanishes (Gross and Kohn 1990). However, if 
one allows more general potentials, the surface integral need not vanish: consider 
fixing an initial state that leads to a certain asymptotic form of no (r ) . Then one can 
always find potentials which increase sufficiently steeply in the asymptotic region 
such that the surface integral does not vanish. In (Xu and Rajagopal 1985; Dhara 
and Ghosh 1987) several examples of this are discussed, where the r.h.s. of Eq. 4.21 
can be zero even while the term inside the divergence is non-zero. These cases are 
however largely unphysical, e.g. leading to an infinite potential energy per particle 
near the initial time. It would be desirable to prove the one-to-one mapping under a 
physical condition, such as finite energy expectation values. 

An interesting case is that of extended periodic systems in a uniform electric 
field, such as is often used to describe the bulk of a solid in a laser field. Repre- 
senting the field by a scalar linear potential is not allowed because a linear potential 
is not an operator on the Hilbert space of periodic functions. The periodic boundary 
conditions may be conveniently modelled by placing the system on a ring. Let us 
first consider a finite ring; for example, a system such as a nanowire with peri- 
odic boundary conditions in one direction, and finite extent in the other dimensions. 
Then an electric field going around the ring cannot be generated by a scalar poten- 
tial: according to Faraday's law, V x E = — dB/dt. Such an electric field can 
only be produced by a time-dependent magnetic field threading the center of the 
ring. The situation for an infinite periodic system in a uniform field also requires 
a vector potential if one wants the Hamiltonian to preserve periodicity. The vector 
potential is purely time-dependent in this case, and leads to a uniform electric field 
via E(t) = — (l/c)dA(t)/dt. Once again, TDDFT cannot be applied, even though 
B = V x A = O.That is, for a finite ring with an electric field around the ring, there 
is a real, physical magnetic field, while for the case of the infinite periodic system 
there is no physical magnetic field, but a vector potential is required for mathe- 
matical reasons. In either case, TDDFT does not apply (Maitra 2003). However, 
fortunately, the theorems of TDDFT have been generalized to include vector poten- 
tials (Ghosh and Deb 1988), leading to time-dependent current-density functional 




(4.32) 
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theory (TDCDFT) (Vignale and Kohn 1996), which will be discussed in detail in 
Chap. 24. Moreover, if the optical response is instead obtained via the limit q —> 0, 
the problem can be formulated as a scalar field and TDDFT does apply. In this case, 
the surface in the RG proof can be chosen as an integral multiple of q and the period- 
icity of the system, and the surface term vanishes. Finally, note that a uniform field 
does not usually pose any problem for any finite system, where the asymptotic decay 
of the initial density in physical cases is fast enough to kill the surface integral. 



4.4.2 Interacting and Non-interacting v-Representability 

The RG proof presented above proves uniqueness of the potential that generates 
a given density from a given initial state, but does not prove its existence. The 
question of whether a given density comes from evolution in a scalar potential is called 
v-representability, a subtle issue that arises also in the ground-state case (Dreizler 
and Gross 1990). In both the ground-state and time-dependent theories, it is still an 
open and difficult one. (See also Fig. 4.1.) Some discussion for the time-dependent 
case can be found in (Kohl and Dreizler 1986; Ghosh and Deb 1988). 

Perhaps more importantly, however, is the question of whether a density, known 
to be generated in an interacting system, can be reproduced in a non-interacting one. 
That is, given a time-dependent external potential and initial state, does a KS system 
exist? This question is called "non-interacting v-representability" and was answered 
under some well-defined conditions in (van Leeuwen 1999). A feature of this proof 
is that it leads to the explicit construction of the KS potential. Chapter 9 covers this 
in detail. 

One condition is that the density is assumed to be time-analytic about the initial 
time. The Runge-Gross proof only requires the potential to be time- analytic, but 
that the density is also is an additional condition required for the non-interacting 
v-representability proof of (van Leeuwen 1999). In fact, it is a much more restrictive 
condition than that on the potential, as has been recently discussed in (Maitra 2010). 
The entanglement of space and time in the time-dependent Schrodinger equation 
means that spatial singularities (such as the Coulomb one) in the potential can lead 
to non-analyticities in time in the wavefunction, and consequently the density, even 
when the potential is time- analytic. Again, we stress that non time-analytic densities 
are covered by the RG proof for the one-to-one density-potential mapping (Sect. 4.2). 
Two different non- time- analytic densities will still differ in their formal time-Taylor 
series at some order at some point in space, and this is all that is needed for the proof. 

The other conditions needed for the KS-existence proof are much less severe. The 
choice of the initial KS wavefunction is simply required to satisfy Eqs. 4.29 and 4.30 
(and be well-behaved in having finite energy expectation values). 

It is interesting to point out here that much less is known about non-interacting 
t>-representability in the ground- state theory. 
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Fig. 4.1 A cartoon to illustrate the RG mapping. The three outer ellipses contain possible density 
evolutions, n(r,t), that arise from evolving the initial state labelling the ellipse in a one-body 
potential t> ext ; the potentials are different points contained in the central ellipse. Different symbols 
label different n(r,t) and one may find the same n(r, t) may live in more than one ellipse. The RG 
theorem says that no two lines from the same ellipse may point to the same v ext in the central ellipse. 
Lines emanating from two different outer ellipses may point to either different or same points in 
the central one. If they come from identical symbols, they must point to different points (i.e. two 
different initial states may give rise to the same density evolution in two different potentials). The 
initial state labelling the lower ellipse has either a different initial density or current than the initial 
states labelling the upper ellipses; hence no symbols inside overlap with those in the upper. However, 
they may be generated from the same potential v ext . Some densities, the open symbols, do not point 
to any u ext , representing the non-u-representable densities. Finally, if an analogous cartoon was 
made for the KS system, whether all the symbols that are solid in the cartoon above, remain solid 
in the KS cartoon, represents the question of non-interacting u-representability 



4.4.3 A Variational Principle 

It is important to realise that the TDKS equations do not follow from a variational 
principle: as presented above, all that was needed was (i) the Runge-Gross proof that a 
given density evolving from a fixed points to a unique potential, for interacting and 
non-interacting systems, and (ii) the assumption of non-interacting v-representability. 

Nevertheless, it is interesting to ask whether a variational principle exists in 
TDDFT. In the ground-state, the minimum energy principle meant that one need only 
approximate the xc energy as a functional of the density, E xc [n], and then take its 
functional derivative to find the ground-state potential: v^[n](r) = 8E xc [n]/8n(r). 
Is there an analogue in the time-dependent case? Usually the action plays the 
role of the energy in time-dependent quantum mechanics, but here the situation 
is not as simple: if one could write v xc [n](r,t) = 8A xc [n]/8n(r, t) for some xc 
action functional A xc [ft] (dropping initial-state dependence for simplicity for this 
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argument), then we see that 8v xc [n](r, t)/8n(r\ t') = 8 2 A xc [n]/8n(r, t)8n(r f , f') 
which is symmetric in t and ^. However that would imply that density-changes at 
later times t' > t would affect the xc potential at earlier times, i.e. causality would 
be violated. This problem was first pointed out in (Gross et al. 1996). Chapter 9 
discusses this problem, as well as its solution, at length. Indeed one can define an 
action consistent with causality, on a Keldysh contour (van Leeuwen 1998), using 
"Liouville space pathways" (Mukamel 2005) or using the usual real-time definition 
but including boundary terms (Vignale 2008). 



4.4.4 The Time-Dependent Current 

Step 1 of the RG proof proved a one-to-one mapping between the external poten- 
tial and the current-density, while the second step invoked continuity, with the help 
of a surface condition, to prove the one-to-one density -potential mapping. In fact, 
a one-to-one mapping between current-densities and vector potentials, a special 
case of which is the class of scalar potentials, has been proven in later work (Xu 
and Rajagopal 1985; Ng 1989; Ghosh and Deb 1988; Gross et al. 1996). But as 
will be discussed in Chap. 24, even when the external potential is merely scalar, 
there can be advantages to the time-dependent current-density functional theory 
(TDCDFT) framework. Simpler functional approximations in terms of the time- 
dependent current-density can be more accurate than those in terms of the density: in 
particular, local functionals of the current-density correspond to non-local density- 
functionals, important in the optical response of solids for example, and for polariz- 
abilities of long-chain polymers. 

In TDCDFT, the KS system is defined to reproduce the exact current-density 
j(r, t) of the interacting system, but in TDDFT the KS current jKs(f% is not 
generally equal to the true current. As both current densities satisfy the continuity 
equation with the same density, i.e. 



we know immediately that the longitudinal parts of j and j'ks must be identical. 
However, they may differ by a rotational component: 



h(r, t) = -V • j(r, t) = -V • jKs(r, t), 



(4.33) 



j(r,t) = j K s(r,t) + j xc (r,t) 



(4.34) 



where 



jxc(r, t) = V x C(r, t) 



(4.35) 



with some real function C(r,t). We also know for sure that 




(4.36) 
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because J d 3 r j(r,t) = J d 3 rrn(r, t) and the KS system exactly reproduces the true 
density. 

The question of when the KS current equals the true current may equivalently 
be posed in terms of v-representability of the current: Given a current gener- 
ated by a scalar potential in an interacting system, is that current non-interacting 
v-representable? That is, does a scalar potential exist in which a non-interacting 
system would reproduce this interacting, v-representable current exactly? By Step 
1 of the RG proof applied to non-interacting systems, if it does exist, it is unique, 
and, since it also reproduces the exact density by the continuity equation, it is iden- 
tical to the KS potential. There are two aspects to the question above. First, the 
initial KS Slater determinant must reproduce the current-density of the true initial 
state. Certainly in the case of initial ground-states, with zero initial current, such 
a Slater determinant can be found (Harriman 1981; Zumbach and Maschke 1983), 
but whether one can be always found for a general initial current is open. Assuming 
an appropriate initial Slater determinant can be found, then we come to the second 
aspect: can we find a scalar potential under which this Slater determinant evolves with 
the same current-density as that of the true system? It was shown in (D' Agosta and 
Vignale 2005) that the answer is, generally, no. (On the other hand, non-interacting A - 
representability, in terms of TDCDFT, has been proven under certain time-analyticity 
requirements on the current-density and the vector potentials (Vignale 2004)). In the 
examples of (D' Agosta and Vignale 2005), even if no external magnetic field is 
applied to the true interacting system, one still needs a magnetic field in a non- 
interacting system for it to reproduce its current. 



4.4.5 Beyond the Taylor-Expansion 

The RG proof was derived for potentials that are time-analytic. It does not apply 
to potentials that turn on, for example, like e _c ^ with C > 0,n > 0, or t p with 
p positive non-integer. Note that the first example is infinitely differentiable, with 
vanishing derivatives at t = while higher-order derivatives of the second type 
diverge as t —> + . This is however only a mild restriction, as most potentials are 
turned on in a time-analytic fashion. Still, it begs the question of whether a proof 
can be formulated that does apply to these more general cases. A hope would be 
that such a proof would lead the way to a proof for non-interacting v-representability 
without the additional, more restrictive, requirement needed in the existing proof that 
the density is time-analytic (Sect. 4.4.2). 

A trivial, but physically relevant, extension is to piecewise analytic potentials, for 
example, turning a shaped laser field on for some time T and then off again. These 
potentials are analytic in each of a finite number of intervals. The potential need not 
have the same Taylor expansion in one interval as it does in another, so the points 
where they join may be points of nonanalyticity. It is straightforward to extend the 
RG proof given above to this case (Maitra et al. 2002b). 
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There are three extensions of RG that go beyond the time-analyticity requirement 
on the potential. The first two are in the linear response regime. In the earliest (Ng 
and Singwi 1987), the short-time density response to "small" but arbitrary potentials 
has been shown to be unique under two assumptions: that the system starts from a 
stationary state (not necessarily the ground state) of the initial Hamiltonian and that 
the corresponding linear density -response function is t- analytic. In the second (van 
Leeuwen 2001), uniqueness of the linear density response, starting from the elec- 
tronic ground-state, was proven for any Laplace-transformable (in time) potential. 
As most physical potentials have finite Laplace transforms, this represents a signifi- 
cant widening of the class of potentials for which a 1:1 mapping can be established 
in the linear-response regime, from an initial ground state. This approach is further 
discussed in Chap. 9, where also a third, completely new way to address this question 
is presented via a global fixed-point proof: the one-to-one density-potential mapping 
is demonstrated for the full non-linear problem, but only for the set of potentials that 
have finite second-order spatial derivatives. 

The difficulty in generalizing the RG proof beyond time-analytic potentials, was 
discussed in Maitra et al. (2010), where the questions of the one-to-one density- 
potential mapping and v-representability were reformulated in terms of uniqueness 
and existence, respectively, of a particular time-dependent non-linear Schrodinger 
equation (NLSE). The particular structure of the NLSE is not one that has been studied 
before, and has, so far, not resulted in a general proof, although Chap. 9 discusses new 
progress in this direction (Ruggenthaler 2011b). On a lattice, the NLSE reverts to a 
system of nonlinear ordinary differential equations; this has been exploited in Tokatly 
(2011b) to prove existence and uniqueness of a TDCDFT for lattice systems. The 
framing of the fundamental theorems of time-dependent density-functional theories 
in terms of well-posedness of a type of NLSE first appeared in Tokatly (2009) where 
it arises naturally in Tokatly 's Lagrangian formulation of TD current-DFT, known as 
TD-deformation functional theory (see Chap. 25). The traditional density-potential 
mapping question is avoided in TD-deformation functional theory, where instead 
this issue is hidden in the existence and uniqueness of a NLSE involving the metric 
tensor defining the co-moving frame. Very recently the relation between the NLSE 
of TDDFT and that of TD-deformation functional theory has been illuminated in 
Tokatly (2011a). 



4.4.6 Exact TDKS Scheme and its Predictivity 

The RG theorem guarantees a rigorous one-to-one correspondence between time- 
dependent densities and time-dependent external potentials. The one-to-one corre- 
spondence holds both for fully interacting systems and for non-interacting particles. 
Hence there are two unique potentials that correspond to a given time-dependent 
density n(r, t): one potential, v QXi [n, #b]( r > t), that yields n(r,t) by propagating the 
interacting TDSE with it with initial state &o , and another potential, i>ks [n, @o](r, t), 
which yields the same density by propagating the non-interacting TDSE with the 
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initial state @o- 

Uexttrc, &o](r, t) n(r, t) l^X v KS [n, &oKr, t). (4.37) 
In terms of these two rigorous mappings, the exact TD xc functional is defined as: 





■ /dV^ 
J \r-r'\ 



v xc [n; #b, ®o](r, t) = v KS [n, ®o](r, t) - v ext [n, V ](r, t) 

(4.38) 

In the past years, the exact xc potential Eq. 4.38 has been evaluated for a few (simple) 
systems (Hessler et al. 2002; Rohringer et al. 2006; Tempel et al. 2009; Helbig et al. 
2009; Thiele et al. 2008; Lein and Kummel 2005). The purpose of this exercise is to 
assess the quality of approximate xc functionals by comparing them with the exact 
xc potential (Eq. 4.38). 

Normally, one has to deal with the following situation: We are given the external 
potential vf^ n (r, t) and the initial many-body state ^o- This information specifies 
the system to be treated. The goal is to calculate the time-dependent density from a 
TDKS propagation. The question arises: can we do this, at least in principle, with 
the exact xc functional, i.e. can we propagate the TDKS equation 



id t (pj(r,t) = 



■«£"W) 



/ 



o ,n(r f ,t) 



<Pj(r,t) (4.39) 



with the exact xc potential given by Eq. 4.38? In particular, what is the initial potential 
with which to start the propagation? 

To answer this, we must understand how the functional-dependences in Eq. 4.38 
look at the initial time. To this end, we differentiate the continuity equation Eq. 4.20 
once with respect to time, and find (van Leeuwen 1999) 

n(r, t) = V • [n(r, t)Vv ext (r, t)] - V ■ a(r, t), (4.40) 

where 

a(r, t) = -i(V(t)\\j(r), f + V^Wit)). (4.41) 



At the initial time, Eq. 4.40 shows that as a functional of the density and initial state, 
Vexdn, &o](r,t = 0) depends on n(r, 0), n(r, 0), and #(0) (through Eq. 4.41). 
Although n(r, 0) is determined by the given ^(0) [and so is h(r, 0) ], ri(r, 0) is not. 
Since, at the initial time, we cannot evaluate the time-derivatives to the left [i.e. via 
n(t),n(t — At),n(t — 2 At)], the start of the propagation may appear problematic. 
But, in fact, we are given the external potential by the problem at hand, as in usual 
time-dependent quantum mechanics. Hence in Eq. 4.39 the initial external potential 
is known, the Hartree potential is specified since the initial wavefunction determines 
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the initial density, and the remaining question is: what is the functional dependence 
of the initial xc potential? If, like v ext , this depends on n, then we would have a 
problem starting the TDKS propagation. Fortunately, it does not. 

Applying Eq. 4.40 to the KS system, we replace v QXt with uks in the first term on 
the right, &(t) with <P(t) in Eq. 4.41 and put u ee = there to obtain: 

n(r, t) = V • [n(r, t)Vv KS (r, t)]-V- a KS (r, t), (4.42) 

where 

a K s(r, t) = -i(0(t)\[j(r), f]\<P(t)). (4.43) 
Subtracting Eq. 4.42 from 4.40 for the interacting system, we obtain 



V • 



n(r,t)W 



' f i ,n(r',t) 
J \r-r f \ 



V-[a KS (r,t)-a(r,t)]. (4.44) 



Evaluating this at t = 0, we see that v xc depends only on the initial states, ^(0) and 
(0) . No second-derivative information is needed, and we can therefore propagate 
forward. 

Analysis of subsequent time-steps was done in Maitra et al. (2008) and showed 
that at the Mi time- step, v xc is determined by the initial states and by densities only at 
previous times as expected from the RG proof. This shows explicitly that the TDKS 
scheme is predictive. 



4.4.7 TDDFT in Other Realms 



A number of extensions of the time-dependent density functional formalism to phys- 
ically different situations have been developed. In particular, for spin-polarized 
systems (Liu and Vosko 1989): one can establish a one-to-one mapping between 
scalar spin-dependent potentials and spin-densities, for an initial non-magnetic 
ground state. Extension to multicomponent systems, such as electrons plus nuclei 
can be found in Li and Tong (1986); Kreibich et al. (2004) and is further discussed 
in Chap. 12, to external vector potentials (Ghosh and Deb 1988; Ng 1989) further 
discussed in Chap. 24, and to open systems, accounting for coupling of the elec- 
tronic system to its environment (Burke 2005c; Yuen-Zhou et al. 2010; Appel and 
Di Ventra 2009; Di Ventra and D' Agosta 2007), covered in Chaps. 10 and 1 1. Other 
extensions include time-dependent ensembles (Li and Li 1985), superconducting 
systems (Wacker 1994), and a relativistic two-component formulation that includes 
spin-orbit coupling (Wang 2005b; Peng et al. 2005; Romaniello and de Boeij 2007). 
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4.5 Frequency-Dependent Linear Response 

In this section we derive an exact expression for the linear density response n\ (r , co) 
of an TV-electron system, initially in its ground-state, in terms of the Kohn-Sham 
density -response and an exchange-correlation kernel. This relation lies at the basis 
of TDDFT calculations of excitations and spectra, for which a variety of efficient 
methods have been developed (see Chap. 7 and also (Marques 2006b)). In fact one of 
these methods follows directly from the formalism already presented: simply perturb 
the system at t = with a weak electric field, and propagate the TDKS equations for 
some time, obtaining the dipole of interest as a function of time. The Fourier transform 
of that function to frequency- space yields precisely the optical absorption spectrum. 
However, a formulation directly in frequency- space is theoretically enlightening 
and practically useful. After deriving the fundamental linear response equation of 
TDDFT in Sect. 4.5.1 , we then derive a matrix formulation of this whose eigenvalues 
and eigenvectors yield the exact excitation energies and oscillator strengths. 



4.5.1 The Density-Density Response Function 

In general response theory, a system of interacting particles begins in its ground-state, 
and at t = a perturbation is switched on. The total potential is given by 



The response of any observable to 8v ext may be expressed as a Taylor series with 
respect to 8v QXt . In particular, for the density, 



Linear response is concerned with the first-order term n\(r,t), while higher-order 
response formalism treats the second, third and higher order terms (see Sect. 4.6). 

Staying with standard response theory, n\ is computed from the density-density 
linear response function x as 



v ex t(r, t) = v ext , (r) + 8v ext (r, t) 



(4.45a) 



8v ext (r, = for t < 0. 



(4.45b) 



n(r, t) = ftGs(f) +fli(r, + ^(r, t) + . . . 



(4.46) 



ni(r,t) = 



S?1 



d 3 r' X (rt,r't')8v ext (r',t'). 



(4.47) 



where 



^ext,0 



(4.48) 
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Ordinary time-dependent perturbation theory in the interaction picture defined with 
respect to u e xt,o yields (Wehrum and Hermeking 1974; Fetter and Walecka 1971) 

X (rt, r't') = -iO(t - t f ){V \[fiH (r, t), h Ho (r\ t')]\V ) (4.49) 

where hu = e 1H ° l n(r)e~ 1Hot and 0(x) = 0(1) for r < (>)0 is the step function. 
The density operator is h(r) = X^=i ~ (Note that the presence of the step- 
function is a reflection of the fact that v QXt [n](r, t) is a causal functional, i.e. the 
potential at time t only depends on the density at earlier times t' < t.) Inserting the 
identity in the form of the completeness of interacting states, ^ 7 | = 1, and 

Fourier-transforming with respect tot — t' yields the "spectral decomposition" (also 
called the "Lehmann representation"): 



X(r,r f ,a)) = ^ 



mHr)\Vi)(Vi\n(r')\V ) (Vo\n(r f )\Vi)(Vi\n(r)m) 



(4.50) 

where the sum goes over all interacting excited states ^ , of energy Ej = Eo + Qi , 
with Eo being the exact ground-state energy of the interacting system. This interacting 
response function x is clearly very hard to calculate so we now turn to TDDFT to 
see how it can be obtained via the noninteracting KS system. 

First, the initial KS ground-state: For t < 0, the system is in its ground-state and 
we take the KS system also to be so. The initial density ^gs(^) can be calculated 
from the self-consistent solution of the ground state KS equations 



— + Vext,o(r) + J d-V p — — + v xc [n GS ](r) 



<pf\r)=e j <pf\r) (4.51) 



and 



n G s(r)= X l^fWI 2 - ( 4 - 52 > 

lowest 

Adopting the standard response formalism within the TDDFT framework, 
we notice several things. Because the system begins in its ground-state, there is 
no initial-state dependence (see Sect. 4.3), and we may write n(r,t) = n[v ex t](r, t). 
Also, the initial potential f e xt,0 is a functional of the ground-state density tigs> so 
the same happens to the response function x = Xt^GsL Since the time-dependent 
KS Eqs. 4.26-4.28 provide a formally exact way of calculating the time-dependent 
density, we can compute the exact density response n\(r,t) as the response of the 
non-interacting KS system: 

m(r, t) = r&t' [ dVxKs(", r't')&VKs(r', t f ), (4.53) 



where 8v^s is the effective time-dependent potential evaluated to first order in the 
perturbing potential, and XKs( r t, r't') is the density-density response function of 
non-interacting particles with unperturbed density ^gs- 
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, , 8n(r, t) 
XKs(rt,r't') = 



SvKs(r',t') 



(4.54) 

^KS [«GS] 



Substituting in the KS orbitals cp^ (calculated from Eq. 4.5 1) into Eq. 4.50 we obtain 
the Lehmann representation of the KS density-response function: 

(0)*/ \ (0), x (0)*. / x (0), 

Z<P k (r)<p) \r)(p) (r f )(p K k \r f ) 
(fk~ fj)Sa k aj- J ~ ' , * , (4.55) 
r,^v kj 3 co-(£j -e k ) + iri 

where fk , fj are the usual Fermi occupation factors and denotes the spin orien- 
tation of the &th orbital. 

The KS density-density response function Eq. 4.55 has poles at the bare KS single- 
particle orbital energy differences; these are not the poles of the true density-density 
response function Eq. 4.50 which are the true excitation frequencies. Likewise, the 
strengths of the poles (the numerators) are directly related to the optical absorption 
intensities (oscillator strengths); those of the KS system are not those of the true 
system. We now show how to obtain the true density -response from the KS system. 
We define a time-dependent xc kernel f xc by the functional derivative of the xc 
potential 



, r v . 8v xc [n](r,t) 

fx C [ncs](rt,r t ) = 



8n(r f , t f ) 



(4.56) 

n=ncs 



evaluated at the initial ground state density ^gs • Then, for a given 8 v ext , the first-order 
change in the TDKS potential is 

f , ,ni(r f ,t) 

8v KS (r, t) = 8v QXt (r, t) + / dV^^-f 



+ J &h'J dV/xctnGsKrr.rY^iCr'.r 7 )- (4.57) 

Equation 4.57 together with Eq. 4.53 constitute an exact representation of the linear 
density response (Petersilka 1996a, b). These equations were postulated (and used) in 
Ando (1977a, b); Zangwill and Sovea (1980a, b); Gross and Kohn (1985) prior to their 
rigorous derivation in Petersilka (1996a). That is, the exact linear density response 
n\(r,t) of an interacting system can be written as the linear density response of 
a non-interacting system to the effective perturbation 8v^s(r, t). Expressing this 
directly in terms of the density-response functions themselves, by substituting Eq. 
4.57 into Eq. 4.53 and setting it equal to Eq. 4.47, we obtain the Dyson-like equation 
for the interacting response function: 
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dhj d 3 nj dt 2 J d 3 r 2 XKsD*GsK^,^i) 



+ Acl>*GS](>^l,^2) 



X[nGs](r2t2,r r t r ). 



\r\ -r 2 \ 

(4.58) 

This equation, although often translated into different forms (see next section), plays 
the central role in TDDFT linear response calculations. 



4.5.2 Excitation Energies and Oscillator Strengths 
from a Matrix Equation 

We now take time-frequency Fourier- transforms of Eqs. 4.56-4.58 to move towards a 
formalism directly in frequency-space. The objective is to set up a framework which 
directly yields excitation energies and oscillator strengths of the true system, but 
extracted from the KS system. We write 

X (co) = Xks (co) + Xks (co) * /hxc (co) * X (co) (4.59) 

where we have introduced a few shorthands: First, we have defined the Hartree-xc 
kernel: 

/Hxc(r, r f , co) = — !— + /*.(#•, r', co). (4.60) 
\r-r'\ 

Note that / X cD?Gs](f% r f , co) is the Fourier-transform of Eq. 4.56; the latter depends 
only on the time difference (t — t f ), like the response functions, due to time- 
translation invariance of the unperturbed system, allowing its frequency-domain 
counterpart to depend only on one frequency variable. Second, we have dropped the 
spatial indices and introduced the shorthand ★ to indicate integrals like xks(&0 * 
fnxcico) = / d 3 nxKs(r,ri,co)fuxc(ri,r f ,co) thinking of x,Xks,/hxc etc as 
infinite-dimensional matrices in r,r f , each element of which is a function of co. 
Now, integrating both sides of Eq. 4.59 against 8v QXt (r, co), we obtain 

1 - Xks(co) * /hxc(g>)] *wi(to) = XksO) + 8v ext (co). (4.61) 

The exact density -response n \ (r , co) has poles at the true excitation energies Q\ . 
However, these are not identical with KS excitation energies s a — Si where the poles 
of xks lie, i.e. the r.h.s. of Eq. 4.61 remains finite for co — > Q\ . Therefore, the integral 
operator acting on n\ on the l.h.s. of Eq. 4.61 cannot be invertible for co — >► as it 
must cancel out a pole in n \ in order to create a finite r.h.s. The true excitation energies 
Qi are therefore precisely those frequencies where the eigenvalues of the integral 
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operator on the left of Eq. 4.61, j^l — xks(&0 * /hxc(<^)J , vanish. Equivalently, the 
eigenvalues k(co) of 

Xks(<») * fmdco) * §(©) = X(co)^(co) (4.62) 

satisfy k(&i) = 1 . This condition rigorously determines the true excitation spectrum 
of the interacting system. 

For the remainder of this subsection, we will focus on the case of spin- saturated 
systems: closed-shell singlet systems and their spin-singlet excitations, to avoid 
carrying around too much notation. We shall return to the more general spin- 
decomposed version of the equations in Sect. 4.5.4. 

Before we continue to cast Eq. 4.62 into a matrix form from which excitation 
energies and oscillator strengths may be conveniently extracted, we mention two very 
useful approximations that can shed light on the workings of the response equation. 
The idea is to expand all quantities in Eq. 4.62 around one particular KS energy 
difference, co q = s a — Si say, keeping only the lowest-order terms in the Laurent 
expansions (Petersilka 1996a, b). This is justified for example, in the limit that the 
KS excitation of interest is energetically far from the others and that the correction 
to the KS excitation is small (Appel et al. 2003). This yields what is known as the 
single-pole approximation (SPA), which, for spin- saturated systems, is: 

Q=co q +2V\tJ d 3 r J d 3 r f <P*(r)f mc (r, r' , a> q )<P q (r') (4.63) 

where we have defined the transition density 

<P q (r) = (p*(r)(p a (r). (4.64) 

This approximation is equivalent to neglecting couplings with all other excitations. 
If also the pole at — co q is kept (the backward transition), i.e. 



XKS 



<P*(r f )<P q (r) <P q (r f )<P*(r)' 
CO — co q + i0+ CO + co q + i0+ 



(4.65) 



(where again the factor of 2 arises from assuming a spin- saturated system), then we 
obtain the "small-matrix approximation"(SMA): 

Q 1 = co\ + 4co q J d 3 r J d 3 r f q (r)f mc (r, r' , Q)<P q (r f ) (4.66) 

provided the KS orbitals are chosen to be real. Discussion on the validity of these 
approximations and their use as tools to analyse full TDDFT spectra can be found 
in Appel et al. (2003). Generalized frequency-dependent, or "dressed" versions of 
these truncations have been used to derive approximations in certain cases, e.g. for 
double-excitations (see Chap. 8). 
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We return to finding a matrix formulation of Eq. 4.62 to yield exact excitations 
and oscillator strengths; we follow the exposition of Grabo et al. (2000). For a single 
KS transition from orbital i to a, we introduce the double-index q = (i,a), with 
transition frequency 

co q =s a - Si (4.67) 

and transition density as in Eq. 4.64. Let a q = fa — f a be the difference in their 
ground- state occupation numbers (e.g. a q = if both orbitals are occupied or if both 
are unoccupied, a q = 2 if i is occupied while a unoccupied (a "forward" transition), 
a q = —2 if a is occupied while i unoccupied (a "backward" transition)). Reinstating 
the spatial dependence explicitly and defining 

S q {(o) = J d 3 r ff ^(r f )f mc (r\r ff ,co^(r ff ,co) (4.68) 

we can recast Eq. 4.62 as 

Za a <t> a (r) 
q q , S q (a>) = HcoMr, co). (4.69) 
co — co a + i0 + 

Solving this equation for i=(r,co), and reinserting the result on the r.h.s. of Eq. 4.68, 
we obtain 

^ M aa '{co) 

$ q ,(co)=X(co)$ q (co) (4.70) 



" CO — CO a ' + i0+ 

q * 



where we have introduced the matrix elements 

M qqf {cD)=a qf J d 3 r J &\'cp*(r)f^c{r,r f ,co)cP ql {r f ). (4.71) 

Defining now fi q = $ q (fi)/(fi — co q ), and using the condition that X(£2) = 1 at a 
true excitation energy, we can write, at the true excitations: 

\M qqt {Q) + co q 8 qq r] fi q , = £2p q . (4.72) 

q' 

This eigenvalue problem rigorously determines the true excitation spectrum of 
the interacting system. The matrix is infinite-dimensional, going over all single- 
excitations of the KS system. In practice, it must be truncated. If only forward tran- 
sitions are kept, this is known as the "Tamm-Dancoff " approximation. 

The first matrix formulation of TDDFT linear response was derived in Casida 
(1995, 1996) by considering the response of the KS density matrix. Commonly 
known as "Casida's equations", these equations are similar in structure to TDHF, 
and are what is coded in most of the electronic structure codes today. By considering 
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the poles and residues of the frequency-dependent polarizability, in Casida (1995) 
it was shown that the true frequencies Q\ and oscillator strengths can be obtained 
from eigenvalues and eigenvectors of the following matrix equation: 

RFi = Q\Fi, (4.73) 

where 

Rq q > = co\& qq < + / & r J dV^(r)/Hxc(r, r', ^/)<Zy (r') (4.74) 

The oscillator strength of transition / in the interacting system, defined as 

// = 2 -Qi (|<^ |*l^/>l 2 + I W>I5^/>I 2 + l(^ol£|^/>| 2 ) (4.75) 
can be obtained from the eigenvectors Fj via 

// = 2 - (ixS- 1 / 2 ^! 2 + |v§- 1/2 F 7 | 2 + |zS- 1/2 F/| 2 ) /m 2 (4.76) 

where S//,*/ = S ifk 8jj/a q o) q > with q = (k, /) here. 

The KS orbitals are chosen to be real in this formulation. Provided that real orbitals 
are also used in the secular equation (4.72), Casida's equations and Eq. 4.72 are easily 
seen to be equivalent. The SPA and SMA approximations, Eqs. 4.63 and 4.66 derived 
before, can be readily seen to result from keeping only the diagonal element in the 
matrix M, or in matrix R: neglecting off-diagonal terms, we immediately obtain 
Eq. 4.66. Assuming additionally that the correction to the bare KS transition is itself 
small, we take a square-root of Eq. 4.66 keeping only the leading correction, and 
find the single-pole- approximation Eq. 4.63. 

The matrix equations are often re- written in the literature as 

(™)(?)"(o- '?)(?) 



where 



Aiajb = Sij8 ab (s a - £i ) + 2 Jd 3 rj dV<£*(r)/Hxc(r, r')<P q ,{r') (4.78a) 
Biajb = 2 J d 3 r J d 3 r / 0*(r)/ H xc(r,r / )^-^(r / ), (4.78b) 



which has the same structure as the eigenvalue problem resulting from time- 
dependent Hartree-Fock theory (Bauernschmitt and Ahlrich 1996a). However, we 
note that this form really only applies when an adiabatic approximation (see Sect. 4.7) 
is made for the kernel (but complex orbitals may be used). 

We note that TDDFT linear response equations can be shown to respect the 
Thomas-Reiche-Kuhn sum-rule: the sum of the oscillator strengths equals the 
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number of electrons in the system (see also Chap. 5). This of course is true for the KS 
oscillator strengths as well as the linear-response-corrected ones. As mentioned 
earlier, in the Tamm-Dancoff approximation, backward transitions are neglected, 
e.g. the B-matrix in Eq. 4.77 is set to zero, with the resulting equations resem- 
bling configuration-interaction singles (CIS) [see (Hirata 1999)]. In certain cases 
the Tamm-Dancoff approximation turns out to be nearly as good as (or sometimes 
"better" than (Casida et al. 2000, Hirata and Head-Gordon 1999)) full TDDFT, but 
it violates the oscillator strength sum-rule. 
To summarize so far: 

(i) The matrix formulations (Eqs. 4.72 and 4.73) are valid only for discrete spectra 
and hence are mostly used for finite systems, while the original Dyson-like 
integral equation, Eq. 4.58, is usually solved when dealing with the continuous 
spectra of extended systems. To obtain the continuous part of the spectra of 
finite systems (e.g. resonance widths and positions), the Sternheimer approach, 
described in Chap. 7 is often used. 

(ii) To apply the TDDFT linear response formalism, there are evidently two ingre- 
dients. First, one has to find the elements of the non-interacting KS density- 
response function, i.e. use Eq. 4.51 to find all occupied and unoccupied KS 
orbitals living in the ground-state KS potential uks,o- An approximation is 
needed there for the ground-state xc potential. Second, one has to apply the 
xc kernel / X c,for which in practice approximations are also needed. The next 
section discusses the kernel in a little more detail. 



4.5.3 The xc Kernel 

The central functional in linear response theory f xc is simpler than that in the full 
theory, because instead of functionally depending on the density and its history as 
well as the initial-states as v xc must, it depends only on the initial ground-state density. 
The kernel can be obtained from the functional derivative, Eq. 4.56, but often a more 
useful expression is to extract it from Eq. 4.58: one can isolate f xc in Eq. 4.58 by 
applying the inverse response functions in the appropriate places, yielding 

f*c[n G sKrt, rV) = x KS W](^, rV) - /"WlCr*, rY) - (4.79) 

\r — r'\ 

where Xks anc * X -1 stand for the kernels of the corresponding inverse integral 
operators. 

Note that the existence of the inverse density-response operators on the set of 
densities specified by Eqs. 4.45a-4.47 follows from Eq. 4.21 in the RG proof: the 
right-hand side of Eq. 4.21 is linear in the difference between the potentials. Conse- 
quently, the difference between n(r,t) and n'(r, t) is non-vanishing already in first 
order of v (r , t) — v' (r , t) . This result ensures the invertibility of linear response oper- 
ators. The frequency -dependent response operators x ( r > f f ,co) and xks(*% r',co), on 
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the other hand, can be non-invertible at isolated frequencies (Mearns and Kohn 1987; 
Gross and Deb 1988). Recently, the numerical difficulties that the vanishing eigen- 
values of xks(/, r\ co) cause for exact-exchange calculations of spectra have been 
highlighted in (Hellgren and von Barth 2009). The non-invertibility means that one 
can find a non-trivial (i.e. non- spatially constant) monochromatic perturbation that 
yields a vanishing density response. This might appear at first sight to be a counter- 
example to the density-potential mapping of the RG theorem, as it means that two 
different perturbations may be found which have the same density evolution (at least 
in linear response). However, this can only happen if the perturbation is truly mono- 
chromatic, having been switched on adiabatically from t = — oo. If we instead think 
of the perturbation being turned on infinitely- slowly from t = 0, it must have an 
essential singularity in time: (e.g. v ~ lim rj ^Q+ o~^/ tJrla)t ), i.e. the potential is not 
time-analytic about t = 0, and so is a priori excluded from consideration by the RG 
theorem. 

Due to causality, f xc (rt,r't') vanishes for t' > t, i. e., f xc is not symmetric 
with respect to an interchange of t with t f . Consequently, f xc (rt, r't') cannot be a 
second functional derivative 8 2 F xc [n]/8n (r , t)8n(r' ', t') (Wloka 1971), and the exact 
VxcW( r i t) cannot be a functional derivative, in contrast to the static case. (See also 
earlier Sect. 4.4.3). 

Known exact properties of the kernel are given in Chap. 5. These include symmetry 
in exchange of r and r' and Kramers-Kronig relations for f xc (r, r r ,co). These rela- 
tions make evident that frequency-dependence goes hand-in-hand with f xc (r,r f , co) 
carrying a non-zero imaginary part. 

The manipulations leading to the Dyson-like equation can be followed also in 
the ground-state Hohenberg-Kohn theory to yield static response equations. The 
frequency-dependent interacting and non-interacting response functions are replaced 
by the interacting and KS response functions to static perturbations, and the kernel 
reduces to the second functional derivative of the xc energy. It follows that 



8 2 E xc [n] 



lim f xc [n GS ](r, r> , co) = fT c [n GS ](r, r f ) = 
co^o dn(r)8n(r f ) 



(4.80) 



The adiabatic approximation for the kernel used in almost all approximations today 
takes f xc [n G s](r, r' , co) = / x S c atic [ftGs](f% r r ). We return to this in Sect. 4.7. 

Finally, we note that here we have only dealt with the linear response to time- 
dependent scalar fields at zero temperature. The corresponding formalism for systems 
at finite temperature in thermal equilibrium was developed in Ng and Singwi (1987) 
and Yang (1988). For the response to arbitrary electromagnetic fields, some early 
developments were made in Ng (1989), and, more recently, in the context of TDCDFT 
of the Vignale-Kohn functional, in van Faassen and de Baij (2004); Ullrich and Burke 
(2004). 
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4.5.4 Spin-Decomposed Equations 



The linear-response formalisms presented above focussed on closed-shell singlet 
systems and their singlet excitations. To describe singlet-triplet splittings, linear 
response based on the spin-TDDFT of Liu and Vosko (1989) must be used. We distin- 
guish two situations: first is when the initial ground- state is closed-shell and non- 
degenerate, and second when the initial state is an open-shell, degenerate, system. 
Here we focus on the first case, where the equations given above are straightforward 
to generalize, e.g. for the response of the spin-a -density, we have 

n 1(r (r, *>) = X / d3r W'( r . r '> ^*W'(r', a>) (4.81) 

a' 

where u e xt,<r is the spin-dependent external potential and Xa,a f ( r ^ r \ is the 
spin-decomposed density-density response function. The fundamental Dyson-like 
equation, Eq. 4.58, remains essentially the same, as do the matrix equations for 
the excitation energies, only spin-decomposed, with the spin-dependent xc kernel 
defined via 



fxc,oo'\n^,nQ±\(rt,r t ) = 



(4.82) 



The exact equations, Eq. 4.72, generalize to: 

X X [ M q°q>o>(n) + ai qa 8 qq >8 aa >] p qW = np qo (4.83) 

or' q' 

with the obvious spin-generalized forms of the terms, e.g. 

M q<rq , <r ,(a>)=a q , a , J d 3 r J d 3 r f &* qa (r)f Uxc ^(r,r\cD)0 q ^(r f ) (4.84) 



wither = f io - f ao and ^ q o = (p* a (r)(p aa (r). 

For spin-unpolarized ground- states, there are only two independent combinations 
of the spin-components of the xc kernel because the two parallel components are 
equal and the two anti-parallel are equal: 

Ac = - ^ fxc,aa' = -(/xc,tt + Ac,t|) (4.85a) 

GO 1 

G X c = - ^ fxc,aa f = -(/xc,tt ~~ Ac,f|)« (4.85b) 



oo' 



The spin-summed kernel, f xc in Eq. 4.85a, is exactly the xc kernel that appeared in the 
previous section. For example, for the simplest approximation, adiabatic local spin- 
density approximation (ALDA) (see more shortly in Sect. 4.7), for spin-unpolarized 
ground- states 
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f^ DA [n](r,r') = 8(r-r') 



. d 2 [ne^ m (n)] 



dn 2 



G^ DA [n](r,r') = S(r-r') a ^ m( ; n l r)) 

n(r) 



n=n(r) 



(4.86a) 
(4.86b) 



where e^ m (n) is the xc energy per electron of a homogeneous electron gas of 
density n, and a£° m is the xc contribution to its spin-stiffness. The latter measures 
the curvature of the xc energy per electron of an electron gas with uniform density 
n and relative spin-polarization m = (n^ — n±)/(n^ + n±), with respect to m, at 
m = 0: a^ m (n) = 8 2 e^ m (n, m)/8m 2 \ m=0 . 

For closed-shell systems, there is no singlet-triplet splitting in the bare KS eigen- 
value spectrum: every KS orbital eigenvalue is degenerate with respect to spin. 
However, the levels spin- split when the xc kernel is applied. This happens even 
at the level of the SPA applied to Eq. 4.83 (Petersilka 1996b; Grabo et al. 2000): one 
finds the two frequencies 



^1,2 = o) q + £Ke {M p ^ p ^ ± Mp^p^}. 



(4.87) 



Using the explicit form of the matrix elements (Eq. 4.84) one finds, dropping the 
spin-index of the KS transition density, the singlet and triplet excitation energies 
within SPA, 



^singlet = Wq+ 2Dle J d 3 r J dV<P*(r) 



1 



77 + fxc(r , r' , co q ) 



<P q (r f ) 
(4.88a) 



^ tri P let = co q + 29te j d 3 r j d 3 r>*(r)G xc (r, r' , oo q )<P q (r f ). (4.88b) 

This result shows that the kernel G xc represents xc effects for the linear response 
of the frequency-dependent magnetization density m(r,co) (Liu and Vosko 1989). 
In this way, the SPA already gives rise to the singlet-triplet splitting in the excita- 
tion spectrum. For unpolarized systems, the weight of the pole in the spin- summed 
susceptibility (both for the Kohn-Sham and the physical systems) at £? tnplet is exactly 
zero, indicating that these are the optically forbidden transitions to triplet states. 



4.5.5 A Case Study: The He Atom 

In this section, we take a break from the formal theory and show how TDDFT linear 
response works on the simplest system of interacting electrons found in nature, the 
helium atom. 

Recall the two ingredients needed for the calculation: (i) the ground- state KS 
potential, out of which the bare KS response is calculated, and (ii) the xc kernel. 



4 Introduction to TDDFT 



83 



Continuum- 

5s - 
4s - 



3p 



2p 



KS 



Singlet 



Triplet 



ALDA 



TDOEP 

x-only 



TDOEP 
SIC 



Exact 



0.90 



0.85 B 



0.80 



0.75 



0.70 



Fig. 4.2 Excitations in the Helium atom (from Petersilka (2000)): bare KS excitations out of the 
exact He ground-state potential {left), followed by TDDFT with the ALDA kernel, exact-exchange 
(TDOEP x-only) kernel, and self-interaction corrected LDA (TDOEP SIC) kernel, and finally the 
exact {right) 



In usual practice, this means two approximate functionals are needed. However 
the helium atom is small enough that the essentially exact ground-state KS poten- 
tial can be calculated in the following way. A highly- accurate wavefunction calcu- 
lation can be performed for the ground-state, from which the density no§{r) = 
2 J d 3 r / |^Gs(^,^ / )| 2 can be extracted. The corresponding KS system consists 
of a doubly-occupied orbital, cpo(r) = V^GS ( r )/2, so that the KS equation 
Eq. 4.51 can easily be inverted to find the corresponding ground-state KS poten- 
tial fKs(f) = V 2 <^o/(2^o) + £o> where so = —I, the exact ionization potential. 

First, we demonstrate the effect of the xc kernel, by utilizing the essentially 
exact ground-state KS potential, obtained by the above procedure beginning with 
a quantum monte carlo calculation for the interacting wavefunction, performed by 
Umrigar and Gonze (1994). The extreme left of Fig. 4.2 shows the bare KS excita- 
tions co q = s a — Si . We notice that these are already very close to the exact spectrum, 
shown on the extreme right, and always lying in between the true singlet and triplet 
energies (Savin et al. 1998). The middle three columns show the correction due to the 
TDDFT xc kernel for which three approximations are shown. The first and simplest 
is the ALDA of Eq. 4.86a and the other two are orbital-dependent approximations 
which will be explained in Sect. 4.7. For now, we simply note that the bare KS 
excitations are good zeroth order approximations to the true excitations, providing 
an average over the singlet and triplet, while the approximate TDDFT corrections 
provide a good approximation to their spin- splitting. 
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Table 4.1 Singlet(s) and triplet (t) excitation energies of the helium atom (from Petersilka (2000)), 
in atomic units 
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The exact results are from the variational calculation of Kono (1984). The second column shows 
the single-particle excitations obtained out of the exact KS potential, while the third and fourth 
columns show those of the approximate EXX and LDASIC potentials. The fifth and sixth columns 
then apply the respective xc kernels to get the TDDFT approximations 



For most molecules of interest however, the exact ground-state KS potential is 
not available. Using LDA or semi-local GGA's can give results to within a few 
tenths of an eV for low-lying excitations. However, for higher excitations, (semi)local 
approximations run into problems because the LDA potential asympototically decays 
exponentially instead of as — 1/r as the exact potential does, so the higher lying 
bound- states become unbound. There is no Rydberg series in LDA/GGA atoms. 
For our simple helium atom, the situation is severe: none of the excitations are 
bound in LDA, and GGA does not improve this unfortunate situation. Use of a 
ground-state xc potential that goes as —1/r at long-range pulls these excitations 
down from the continuum into the bound spectrum, and, as Table 4.1 shows, can 
be quite accurate. The table shows results using the exact-exchange approximation 
(EXX), and the self-interaction-corrected local density approximation (LDASIC); 
both bare KS excitations as well as the TDDFT values (i.e. corrected by the kernel) 
are shown. These approximations are discussed in detail in Sect. 4.7; to note for 
now, is that the ground- state KS potential in all cases has the correct long-range 
behavior. Notice also that the bare KS excitations are quite accurate; applying the 
kernel (second step) provides a small correction. 
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Fig. 4.3 The exact and LDA KS potentials for the He atom 



Consider Fig. 4.3 which plots the true and LDA xc potentials for the case of 
the helium atom; similar pictures hold for any atom. In the shaded valence region, 
we notice that the LDA xc potential differs from the exact xc potential by nearly 
a constant. This effect is related to the derivative discontinuity, and it was argued 
in Perdew (1985) that this constant has a value (7 + A)/2 where / is the ioniza- 
tion potential and A the electron affinity. The fact that the LDA xc potential runs 
almost parallel to the exact in this region, means that the valence orbitals are well- 
approximated in LDA, while their orbital energies are almost uniformly shifted up. 
This is why excitations, starting from the zeroth-order KS orbital energy differences 
in the valence region are generally approximated well in LDA. However, the rapid 
decay of v]^ A to the zero-asymptote means that the higher-excitation energies are 
underestimated and eventually get squeezed into the continuum. (Unfortunately for 
the case of the He atom, this happens to even the lowest excitations.) 

The top panel of Fig. 4.4 illustrates two effects of the too rapid decay of the 
LDA potential on the optical spectrum: (i) it pushes the valence levels up, so that 
the ionization potential is too low; the onset of the LDA continuum is red-shifted 
compared to the exact, and (ii) there is no Rydberg series in the LDA spectrum; 
instead their oscillator strengths appear in the LDA continuum, but in fact are not 
badly approximated (Wasserman et al. 2003). The reason for this accuracy is due to 
the LDA and true xc potentials running nearly parallel in the valence region: the LDA 
HOMO orbital, out of which the transitions are computed, very well- approximates 
the true HOMO, while the LDA continuum state at energy E = co + / LDA follows 
very closely the exact continuum state at energy E = co + / exact until a distance 
large enough away from the nucleus that the integrand does not contribute due to the 
decay of the HOMO. Noting that KS spectra are not true spectra, the lower panel 
shows the TDDFT-corrected spectrum using ALDA for the xc kernel; although not 
resolving the discrete part of the spectra, the overall oscillator strength envelope is 
not bad. For a detailed discussion, we refer the reader to Wasserman et al. (2003). 



86 



E. K. U. Gross and N. T. Maitra 



Fig. 4.4 Optical absorption 
spectrum in the He atom 
(from Wasserman et al. 
(2003)). The top panel shows 
the bare exact KS and LDA 
KS spectra. The lower panel 
shows the TDDFT ALDA 
spectra {dashed line, from 
Stener et al. (2001)), the 
exact calculations from Kono 
and Hattori (1984) and 
experimental results from 
Samson et al. (1994) 
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The purpose of this case study was to illustrate the workings of TDDFT response 
and investigate only the simplest functional on the simplest system. We note that 
most molecules have many more lower-lying excited states so that GGA's can do 
a much better job for more excitations. We return to the question of functional 
approximations in later sections throughout this book. 

4.6 Higher-Order Response 

Often the simplest way to calculate the non-linear response of a system to an 
external perturbation is via time-propagation. But, like in the linear-response case, 
it can be instructive to perform the non-linear response calculation directly in 
frequency-space. The higher-order terms in Eq. 4.46 can be expressed in terms of 
higher-order density-density response functions: 

n 2 (x) = \\ J &c' J dx"x (2 \x,x',x")Sv(x')Sv(x") (4.89a) 

m(x) = — f dx f f dx" ( dx fff y^ 2 \x, x\ x ff , x fff )8v(x)8v(x ff )8v(x fff ) 

3!7 J J A (4.89b) 

where we used the short-hand x = (r, t) and J dx = J d 3 rf dt. For quadratic 
response, the analogues of Eqs. 4.48 and 4.49 are 
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: m (x,x',x") 



8 2 n(x) 
8v ext (x')8v ext (x") 

-^0(t -t f )0(t 



t ,f )(Vo\[[n H (x), n H (x% n H (x f, )Wo) 



(4.90) 

where the sum goes over all permutations of x , x f , x" . Clearly, the interacting higher- 
order response functions are very difficult to calculate directly and instead we look to 
extract them from the KS response functions and xc kernels. Manipulations similar 
to those in the linear response case, but more complicated, lead to (Gross et al. 1996): 



xVHx 9 xW')=JdyJdy' X jg(x, 



y,y) 



<^Ks(y) 



8v(x f ) 



8vKs(y r ) 



+ j dvxKs(^,v) j dy f J dy ff k xc [n GS ](y,y f ,y ff )x(y f ,x f ) X (y ff ,x ff ) 
+ J dyxKs(x,y) J dy'fti xc [n G s](y,y')x (2) (y',x',x"). 

(4.91) 

Here Xk* = 8 2 n(x) /8v^s (x f )8v^s (x ff ) I is the KS second-order density -response 
function, and 

8 2 v xc [n](r,t) 



k xc [n G s](rt,r f t f ,r"t") = 



8n(r f ,t f )8n(r",t") 
is the dynamical second-order xc kernel. In the adiabatic approximation, 



(4.92) 

n=nos 



8n(r)8n(r f )8n(r ff ) 



C a [«](r, r', r") = , /V,,. , ,„ (4.93) 



with E xc [n] a ground-state xc energy functional. Making Fourier-transforms with 
respect to t — t' and t — t" , we arrive at the Dyson equation 

n 2 (r, co) =- J dto r J d 3 r\d 3 r 2 x^ 2 \r, r l> r 2, to,to — to ! )8v{r\, co)8v(r 2 , to — to r ) 

= X - j dto J d 3 r x d 3 r 2 ^x^s^ r ^ r ^^ r 2^^-^) 8v KS(ri,to)8v K s(r2,oj-co f ) 
+ J d 3 r 3 XKS(r,ri,co)k xc (ri,r 2 ,r 3 ,co,co- to')n\{r 2 , to')ni(r 3 , co - to f )^j 
+ J d 3 nd 3 r 2 XKS(r, ri, co)fn xc (ri, r 2 , to)n 2 (r 2 , to) 

(4.94) 

Likewise, one may work out Dyson-like response equations for the higher-order 
response functions, each time introducing a new higher-order xc kernel. These deter- 
mine the frequency-dependent non-linear response. Sum-over- states expressions for 
the non-interacting KS density -response functions up to third-order may be found 
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in Senatore and Subbaswamy (1987). We also point to Chap. 7, where higher-order 
response is discussed within a Sternheimer scheme. 

Gross et al. (1996) pointed out a very interesting hierarchical structure that the 
TDDFT response equations have. At any order i, 

Hi (co) = Mi (co) + xks(o>) * fttxc(o>) * m (co) (4.95) 

where Mi depends on lower-order density-response (and response-functions up to 
zth order). The last term on the right of Eq. 4.95 has the same structure for all orders. 
If we define the operator 

L(co) = 1 - XKS(<») * MM (4.96) 

then 

L(cD)*m(cD) = Mi(cD) (4.97) 

so L(co) plays a significant role in determining what new poles are generated from 
electron-interaction effects in all orders of response (Elliott 2011). 



4.7 Approximate Functionals 

As noted earlier, the xc potential is a functional of the density, the initial true state, 
and the initial KS state. The exact functional has "memory", that is, it depends on 
the history of the density as well as these two initial states. This is discussed to some 
extent in Chap. 8. In fact these two sources of memory are intimately related, and often 
the elusive initial-state dependence can be replaced by a type of history-dependence. 
The xc kernel of linear response has simpler functional dependence, as it measures 
xc effects around the initial ground- state only. Functionally it depends only on the 
initial ground- state density, while memory-dependence appears as dependence on 
frequency in the arguments of f xc . 

It should be noted that when running response calculations, for formal consistency, 
the same approximation should be used for the xc kernel as is used for the ground- 
state potential, i.e. there must exist an approximate functional v^ v [n](r, t) such that 
the initial potential v^ v [n](r, t = 0) is used in the KS ground-state calculation, 
and such that fxc V [n](r, r f ,t — t') = 8v^ v [n](r, t)/8n(r f , t f ) is used in the time- 
dependent response part. If different functionals were used for each step, one has 
left the framework of TDDFT response, since the calculation no longer is equivalent 
to computing the time-dependent response to an external perturbation. Nevertheless, 
this fact is often ignored in practical calculations. 

We now will outline some of the different approximations people use today. 
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4.7.1 Adiabatic Approximations: ALDA, AGGA, AB3LYP, etc. 

Almost all functionals in use today have no memory-dependence whatsoever. They 
are "adiabatic", meaning that the density at time t is plugged into a ground- state 
functional, i.e. 



,,adia 



[n](r,t) = vg s [n(t)](r). 



(4.98) 



For the xc kernel, we retrieve the static kernel of Eq. 4.80 



A a c dia ["GsK^,rV) 



8vg s [n(t)](r) 



8n(r f , t f ) 



8(t -t f ) 



8 2 E xc [n] 



8n(r)8n(r f ) 



(4.99) 



If the external time-dependence is very slow (adiabatic) and the system begins in 
a ground-state, this approximation is justified. But this is not the usual case. Even 
if the density is reproducible by a system in its a ground-state, the wavefunction is 
usually not, so this appears to be quite a severe approximation. Nevertheless adiabatic 
approximations are the workhorse of the myriads of applications of TDDFT today, 
and work pretty well for most cases (but not all). Why this is so is still somewhat 
of an open question. Certainly adiabatic approximations trivially satisfy many exact 
conditions related to memory-dependence, so perhaps this is one reason. This is 
similar to the justification used for the success of LDA in the ground-state case, and 
the subsequent development of generalized gradient approximations (GGA) based on 
satisfaction of exact conditions. When considering excitation energies of systems, the 
bare KS orbital energy differences themselves are reasonably good approximations 
to the exact excitation energies. The kernel then just adds a small correction on top of 
this good zeroth order estimate, and hence even the simplest approximation such as 
an adiabatic one does a decent job. Many cases where the usual approximations fail, 
such as excited states of multiple-excitation character, or certain types of electronic 
quantum control problems, can be clearly understood to arise from lack of memory 
in the adiabatic approximation. 

The adiabatic local density approximation, or ALDA is the simplest possible 
approximation in TDDFT. It is also often called TDLDA (for time-dependent LDA): 



i4 LDA M(r, t) = v^ m (^(r, 0) = ^[ne^ m (n)]\ n=n(rJ) (4.100) 

an 

where e^ m (n)) is the xc energy per particle of the homogeneous electron gas. The 
corresponding xc kernel 

d 2 

A^ LDA ["gsK^, rV) = 8(t - t f )8(r - r f )^ [ne*° c m (n)]\ n=nGs{r) (4.101) 

which is completely local in both space and time, and its Fourier- transform, Eq. 4.86a 
is frequency-independent. Although it might appear justified only for slowly- varying 
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systems in space and in time, it often gives reasonable results for systems far from 
this limit. Adiabatic GGA's and hybrid functionals are most commonly used for finite 
systems; hybrids in particular for the higher-lying excitations where it is important 
to catch the tail of the molecular potential. 



4 J. 2 Orbital Functionals 

A natural way to break free of the difficulties in approximating functionals of the 
density alone, and still stay within TDDFT, is to develop functionals of the KS 
orbitals. The simplest functional of this kind is the exact-exchange functional, derived 
from the action: 



A x [{0/ a }] 



oo 

3 f cp* a (r\ t')<p jo {r', t')<fH a (r, t f )cp* a (r, *') 



a iJ 
X / dV 



Note that for more general functionals, the action needs to be defined on the Keldysh 
contour (van Leeuwen 1998) (see also Chap. 6). The exact exchange potential is then 
given by 

VxAWjamr, t) = 8 f M ^ ]] (4.103) 

Orbital functionals are in fact implicit density functionals because orbitals are triv- 
ially functionals of the single-particle KS potential, which, by the RG theorem, is a 
functional of the density, cpj [vks] M ir,t). The xc potential is given by the functional 
derivative of the action with respect to the (spin-) density and the xc kernel is the 
second functional derivative. The equation satisfied by the xc potential is usually 
called the (time-dependent) Optimized Effective Potential (OEP) equation, and is 
discussed in Chap. 6. The exact-exchange functional is local in time when viewed as 
a functional of KS orbitals. However, viewed as an implicit functional of the density, 
it is non-local in time. The second- functional derivative with respect to the density 
then has a non- trivial dependence on t — t' . 

Invoking a Slater-type approximation in each functional derivative of Eq. 4.102, 
(Petersilka 1996a, 1998) deduced, 

r') = 1 ^f^f^\\ (4 . 10 4) 

|r-r'| n a (r)n a (r f ) 

Evidently, with this approximation, the non-trivial (t — t') -dependence of the exact 
exchange-only kernel is not accounted for. However, it is clearly spatially non-local. 
For one and two electrons, Eq. 4.104 is exact for exchange. 
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The full numerical treatment of exact-exchange (TDEXX), including memory- 
dependence, has recently seen some progress. In Gorling (1997b), the exact x kernel 
was derived from perturbation theory along the time-dependent adiabatic connection. 
One scales the electron-electron interaction by X, defining a Hamiltonian 

H k = f + AVee + V k (t) (4.105) 

such that the density n k (r,t) = n(r,t) for any X much like as is done in ground- 
state DFT. The initial state is chosen to reproduce the same initial density and 

its first time-derivative, at any X. For X = 0, V k (t) = Vks(0 an d for A = 1, 
V k (t) = Vext (0 • Performing perturbation theory to first order in X yields the TDEXX 
potential and kernel, found in Gorling (1997, 1998a, 1998b) while higher orders give 
correlation functionals. Until very recently there was only limited use of this kernel 
due to numerical instabilities (Shigeta et al. 2006). A series of papers reformulated 
the problem in terms of response of the KS potential itself, avoiding the calcula- 
tion of numerically prohibitive inverse response functions (Hesselmann et al. 2009; 
Gorling et al. 2010; Ipatov 2010), but needing a time-consuming frequency-iteration 
for each excitation energy. Most recently, a very efficient method has been derived 
that translates the problem onto a generalized eigenvalue problem (Hesselmann and 
Gorling 2011). Although the results of full exact-exchange calculations for excita- 
tion energies are often numerically close to those of time-dependent Hartree-Fock 
(Hesselmann and Gorling 2011) in the cases so far studied, there is a fundamental 
difference in the two methods: TDEXX operates with a local (multiplicative) poten- 
tial, while that of time-dependent Hartree-Fock is non-local, i.e. an integral oper- 
ator. Furthermore, the Hartree-Fock single-particle energy differences s^ F — sf F 
are usually too large and hence the Hartree-Fock kernel reduces the Hartree-Fock 
energy difference, while s^ s —sf s within EXX tend to be too small, so the x-kernel of 
TDDFT has to increase the KS excitation energy. Beyond the linear response regime, 
(Wijewardane and Ullrich 2008) computed nonlinear dynamics in semiconductor 
wells within TDEXX. 

Another class of orbital-dependent functionals are self-interaction-corrected (SIC) 
functionals. An approximation at a similar level to Eq. 4.104 can be found in 
Petersilka (2000). 



4.73 Hydrodynamically Based Kernels 

The first proposal to incorporate memory-dependence was that of Gross and Kohn 
(1985), who suggested to use the frequency-dependent xc kernel of the homogeneous 
electron gas in the sense of an LDA: 

A L c D V G s](r, r\ co) = f^ m (n GS (r), \r - r'|, co). (4.106) 

and furthermore that the response n \ (r , co) is slowly varying enough on the length- 
scale of /xc° m (^GS ( r ) > \r —r f \\co) that only its uniform component contributes. That 
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is, taking a spatial Fourier-transform with respect to r — r f , we include only the 
zeroth-Fourier component. This gives the Gross-Kohn kernel: 

f£ K [n G s](r, r', a>) = S(r - r f )f^ m (n GS (r), q=0,o>) (4.107) 

where q is the spatial Fourier transform variable. One requires the knowledge of 
the frequency-dependent response of a uniform electron gas, about which, indeed 
many exact properties are known, and parametrizations, believed to be accurate, exist 
(Conti 1999, 1997; Gross and Kohn 1985; Qian and Vignale- 2002, 2003). Chapter 24 
discusses some of these. 

Although the GK approximation has memory, it is completely local in space, a 
property which turns out to violate exact conditions, such as the zero-force rule and 
translational invariance (see Chaps. 5 and 24). Even ALDA does not violate these. 
To go beyond the adiabatic approximation consistently, both spatial and temporal 
non-locality must be included. This is perhaps not surprising in view of the fact that 
the density that at time t is at location r was at an earlier time t' < t at a different 
location, i.e. memory is carried along with the fluid element. The development of 
memory-dependent functionals, often based on hydrodynamic schemes, is discussed 
further in Chaps. 8 and 24. These include the Dobson-Biinner-Gross (Dobson et al. 
1997; Vignale and Kohn 1996; Tokatly 2005a, b), and (Kurzweil and Baer 2004) 
approaches. These are not commonly used; only the Vignale-Kohn functional has 
seen a few applications. 



4.8 General Performance and Challenges 

As shown in Sect. 4.2, TDDFT is an exact reformulation of non-relativistic time- 
dependent quantum mechanics. In principle, it yields exact electronic dynamics 
and spectra. In practice, its accuracy is limited by the functional approximations 
used. The simplest and computationally most efficient functional, ALDA, is local 
in both space and time, and it is perhaps surprising that it works as well as it does. 
We now discuss cases where it is essential to go beyond this simple approximation, 
and, further, beyond its adiabatic cousins. We organize this section into three parts: 
linear response in extended systems, linear response in finite systems, and real-time 
dynamics beyond the perturbative regime. 



4.8.1 Extended Systems 

We first ask, how well does ALDA perform for the response of solids? In simple 
metals, ALDA does well, and captures accurately the plasmon dispersion curves 
(Quong and Eguilez 1993). In fact, the ordinary plasmon is captured reasonably 
even by the Hartree potential alone, i.e. setting f xc = 0, which is called the RPA. 
Applying f^ LDA improves the description of its dispersion and linewidth. 
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The answer to the question above is rather more subtle for non-metallic systems. 
ALDA does a good job for electron energy loss (EEL) spectra, both when the 
impinging electron transfers finite momentum q and in the case of vanishing 
momentum-transfer. The EEL spectrum measures the imaginary part of the inverse 
dielectric function (see Chap. 3), which, in terms of the density-response function, is 
given in Eq. 3.56a (Onida et al. 2002; Botti et al. 2007). The optical response, which 
measures the imaginary part of the macroscopic dielectric function (see Chap. 3) is 
also quite well predicted by ALDA (Weissker et al. 2006) for finite wavevector q. 

However, for optical response in the limit of vanishing wavevector, q —> 0, 
ALDA performs poorly for non-metallic systems. There are two main problems: 
(i) The onset of continuous absorption is typically underestimated, sometimes by as 
much as 30-50%. This problem is due to the fact that the KS gap in LDA is much 
smaller than the true gap. But even with the exact ground-state potential, there is very 
strong evidence (Knorr and Godlay 1992; Griming et al. 2006, Niquet and Gonze 
2004) that the exact KS gap is typically smaller than the true gap. To open the gap, 
the xc kernel must have an imaginary part (Giuliani and Vignale 2005). This follows 
from the fundamental Dyson equation (4.59) 2 : We know that the imaginary part of 
Xks ( r > r ' > co) = for co inside the KS gap. Then, for an approximate f xc (r,r f , co) that 
is real, taking the imaginary part of Eq. 4.59 for co inside the KS gap (0 < co < E^ s ), 
yields 

3m xks (co) = — ► [l - £K*xks(*>) * fm c (co)] * 3m X (co) = 0. (4.108) 

Following the analogous procedure for co inside the true gap (0 < co < E g ) , where 
3mx(r, r'co) = 0, yields 

3m X (co) = — ► 3m X Ks(co) * [l + /hkcM^xM] = 0. (4.109) 

In view of the fact that fy xc = Xks ~~ X ~\ me expressions inside the square brackets 
in (4.108) and (4.109) cannot vanish identically in the full interval < co < E^ s 
and < co < Eg, respectively. (The expressions may vanish at isolated frequencies 
corresponding to collective excitations). Hence we must conclude that wherever 
3m Xks (co) = 0, then also 3mx(co) = 0. That is, for frequencies inside the KS 
gap, the true response is also zero. Likewise, wherever 3mx(co) = 0, then also 
3m xks (co) = 0. That is, for frequencies inside the true gap, the KS response is also 
zero. Putting the two together implies that the KS system and the true system must 
have the same gap when an approximation for fy xc is used that is purely real. This 
is clearly a contradiction, implying that /hxc must have a non- vanishing imaginary 
part. This, on the other hand, is equivalent to f xc having a frequency-dependence, as 
mentioned in Sect. 4.5.3. Any adiabatic approximation however takes a kernel that 
is the second density-functional-derivative of a ground- state energy functional, and 
therefore is purely real. Further, as we shall shortly discuss, the kernel must have have 



2 This argument is due largely to Giovanni Vignale. 
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a long-ranged part, that goes as l/q 2 as q —> 0, to get any non- vanishing correction 
on the gap. ALDA, on the other hand, is local in space, and so constant in q. 

(ii) The second main problem is that ALDA cannot yield any excitonic structure; 
again, one needs a long-ranged l/q 2 part in the kernel to get any significant improve- 
ment on RPA for optical response. To reproduce excitons, an imaginary part to f xc 
is not required (Reining 2002). 

Fundamentally, the long-ranged behaviour of f xc can be deduced from exact 
conditions satisfied by the xc kernel, such as the zero-force theorem, that inextricably 
link time-nonlocality and space-non-locality. This is shown explicitly in Chap. 24. 
Lack of a long -ranged term in ALDA or AGGA for finite systems, or for EELs 
spectra, is not as critical as for the case of extended systems, since its contribution 
is much smaller there. 

The need for this long-ranged behavior in the xc kernel is, interestingly, not a 
consequence of the long-rangedness of the Coulomb interaction. A simple way to 
see this is to consider the SPA Eq. 4.63 for a system of size L 3 , where, for the 
extended system we consider L —> oo. The transition densities <P q scale as 1/L 3 , 
so for a finite-ranged xc kernel, the xc-correction to the RPA value scales as 1/L 3 , 
and so vanishes in the extended-system limit (Giuliani and Vignale 2005). 

An alternate way of seeing the need for the long-ranged kernel, is to note that 
the optical absorption measures the imaginary part of the macroscopic dielectric 
function, which can be written in terms of a modified density-response function 
(Botti et al. 2007) (and see Eqs. 3.51 and 3.55 in Chap. 3). Now XKs(q -> 0) ~ q 2 
for infinite systems, so if f xc is to have any significant non- vanishing effect on the 
optical response, it must have a component that diverges as l/q 2 as q —> 0. 

Recent years have seen a tremendous effort to confront the problem of optical 
response in solids by including spatially-non-local dependence. Exact-exchange 
was shown in Kim and Gorling (2002) to have apparent success in capturing the 
exciton. However, it was shown later in Bruneval (2006) that if done carefully, 
the excitonic structure predicted by exact-exchange is far too strong, essentially 
collapsing the entire spectrum onto the exciton. The earlier calculation of Kim and 
Gorling (2002) fortuitously induced an effective screening of the interaction, since 
the long wavelength contributions were cut off in those calculations. TDCDFT has 
also been used (de Boeij 2001), with the motivation that local functionals of the 
current-density contain non-local information of the density, and this is discussed 
further in Chap. 24. Perhaps the most intense progress has been made in the devel- 
opment of kernels derived from many -body perturbation theory (MBPT), leading 
to what is now known as the "nanoquanta kernel". The latter is deduced from the 
Bethe-Salpeter approach of MBPT (Bruneval et al. 2005; Reining et al. 2002; Sottile 
et al. 2003; Adragna et al. 2003; Marini et al. 2003b; Stubner et al. 2004; von Barth 
et al. 2005). An important aspect is that the reference systems in the Bethe-Salpeter 
approach and the TDDFT approach are completely different: KS excitation energies 
and orbitals of the latter are not quasiparticle energies and wavefunctions that the 
former builds on. From the point of view of MBPT, the TDDFT xc kernel may be 
interpreted as having two roles: shifting the KS excitations to the quasiparticle ones 
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(so-called / xc ,) and then accounting for the electron-hole interaction (so-called 

/£>)• 

(2) 

In practical uses of the nanoquanta kernel, explicit models for / x y are used on 
top of simply the quasiparticle energies (usually in GW approximation). We refer 
the reader to Botti et al. (2007) and Onida et al. (2002). 



4.8.2 Finite Systems 

In linear response calculations of optical spectra, use of local or semi-local 
functionals for low-lying excitations, or hybrid functionals for higher-lying ones, 
within the adiabatic approximation for the xc kernel yields results that are typi- 
cally considerably better than those from TDHF or configuration interaction singles 
(CIS). These methods scale comparably to TDDFT, while the accuracy of TDDFT is, 
in most cases, far superior. 

Most quantum chemistry applications use the B3LYP hybrid functional (Becke 
1993a, b). While excitation energies are typically good to within a few tenths of an 
eV, structural properties fare much better (Furche 2002a; Elliott et al. 2009). For 
example, bond-lengths of excited states are within 1%, dipole moments and vibra- 
tional frequencies to within 5%. Chapter 16 discusses this more. Often the level of 
accuracy has been particularly useful in explaining, for the first time, mechanisms of 
processes in biologically and chemically relevant systems, e.g. the dual fluorescence 
in dimethyl-amino-benzo-nitrile (Rappoport and Furche 2004), and chiral identifi- 
cation of fullerenes (Furche 2002b). 

In the following, we discuss several cases where the simplest approximations like 
ALDA and AGGA, perform poorly. 

To be able to describe Rydberg excitations, it is essential that the ground-state 
potential out of which the bare KS excitations are computed has the correct — 1/r 
asymptotics. LDA and GGA do not have this feature, and as we have seen already in 
our case study of the He atom, the Rydberg excitations were absent in LDA/GGA. 
Solutions include exact-exchange methods, self-interaction corrected functionals, 
and hybrids. Step-like features in the ground-state potential as well as spatial non- 
locality in the xc kernel can also be essential: a well-known case is in the compu- 
tation of polarizabilities of long-chain molecules (van Faassen 2002; van Gisbergen 
1999b; Gritsenko 2000), and exact-exchange, as well as TDCDFT-methods have 
been explored for this problem. A more challenging problem is that of molecular 
dissociation: it is notoriously difficult to obtain accurate ground-state dissociation 
curves, since self-interaction errors in the usual functionals leads to fractional charges 
at large separation. As it dissociates, the exact ground-state potential for a molecule 
composed of open-shell fragments such as LiH, develops step and peak features in the 
bond-midpoint region (Perdew 1985; Gritsenko 1996; Helbig et al. 2009; Tempel et 
al. 2009), missed in GGAs and hybrids alike, but crucial for a correct description. To 
get even qualitatively correct excited state surfaces, frequency-dependence is crucial 
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in the xc kernel (Maitra 2005b; Maitra and Tempel 2006a) (see also Chap. 8). The 
essential problem is that the true wavefunction has wandered far from a single- Slater- 
determinant, making the work of the ground- state exchange-correlation potential and 
the kernel very difficult. Such cases of strong correlation are one of the major moti- 
vators of time-dependent density-matrix functional theory, discussed in Chap. 26. 

Another case where frequency-dependence is essential are states of double- 
excitation character. We will defer a discussion of this to Chap. 8. 

A notorious failure of the usual approximations for finite systems is for charge- 
transfer excitations at large-separation (Dreuw 2003, 2004; Tozer 2003). To leading 
order in 1 /R , the exact answer for the lowest charge-transfer excitation frequency is: 

^act _^ jD _ A A_ y R (4 nQ) 

where I D is the ionization energy of the donor, A A is the electron affinity of the 
acceptor and — l/R is the first electrostatic correction between the now charged 
species. Charge-transfer excitations calculated by TDDFT with the usual approxi- 
mations however severely underestimate Eq. 4.110. Due to the exponentially small 
overlap between orbitals on the donor and acceptor, located at different ends of 
the molecule, f xc must diverge exponentially with their separation in order to 
give any correction to the bare KS orbital energy difference (see e.g. Eq. 4.63). 
Semilocal functional approximations for f xc give no correction, so their predic- 
tion for charge-transfer excitations reduces to the KS orbital energy difference, 
s a — Si = si (acceptor) — sh (donor), where L, H subscripts indicate the KS LUMO 
and HOMO, respectively. This is a severe underestimate to the true energy, because 
the ionization potential is typically underestimated by the HOMO of the donor, due 
to the lack of the — 1/r asymptotics in approximate functionals (Sect. 4.5.5), while 
the LUMO of the acceptor lacks the discontinuity contribution to the affinity. The 
last few years have seen many methods to correct the underestimation of CT excita- 
tions, e.g. (Autschbach 2009; Tawada et al. 2004; Vydrov 2006; Zhao and Truhlar 
2006; Stein et al. 2009a; Hesselmann et al. 2009; Rohrdanz 2009); most modify 
the ground- state functional to correct the approximate KS HOMO's underestimation 
of / using range- separated hybrids that effectively mix in some degree of Hartree- 
Fock, and most, but not all (Stein et al. 2009a; Hesselmann et al. 2009) determine 
this mixing via at least one empirical parameter. Fundamentally, staying within pure 
DFT, both the discontinuity contribution to A and the — l/R tail in Eq. 4.1 10 come 
from /hxo which must exponentially diverge with fragment separation (Gritsenko 
and Baerends 2004). Worse still, in the case of open-shell fragments, not covered by 
most of the recent fixes, additionally the exact f xc is strongly frequency-dependent 
(Maitra 2005b). 

We briefly mention two other challenges, which are further discussed later in 
this book. The difficulty that usual functionals have in capturing potential energy 
surfaces near a conical intersection is discussed in Chap. 14. This poses a challenge 
for coupled electron-dynamics using TDDFT, given that conical intersections are 
a critical feature on the potential energy landscape, funneling nuclear wavepackets 



4 Introduction to TDDFT 



97 



between surfaces. The second challenge is the Coulomb blockade phenomenon in 
calculations of molecular transport. The critical need for functionals with a deriv- 
ative discontinuity to describe this effect, and to obtain accurate conductances in 
nanostructures, is further discussed in Chap. 17. 



4.8.3 Non-perturbative Electron Dynamics 

In Chap. 18 of this book, we shall come back to the fascinating world of strong-field 
phenomena, several of which the usual approximations of TDDFT have had success 
in describing, and some of which the usual approximations do not capture well. 
Developments in attosecond laser science have opened up the possibility of elec- 
tronic quantum control; recently the equations for quantum optimal control theory 
within the TDKS framework have been established and this is described in Chap. 13. 
Chaps. 14 and 15 discuss the difficult but extremely important question of coupling 
electrons described via TDDFT to nuclear motion described classically, in schemes 
such as Ehrenfest dynamics and surface-hopping. Here instead we discuss in general 
terms the challenges approximations in TDDFT face for real-time dynamics. 

However first, we show how useful a density-functional picture of electron 
dynamics can be for a wide range of processes and questions, via the time-dependent 
electron localization function (TDELF). With the advent of attosecond lasers, comes 
the possibility of probing detailed mechanisms of electronic excitations and dynamics 
in a given process. For example, in chemical reactions, can we obtain a picture of 
bond-breaking and bond-forming? In Burnus (2005) it was shown how to generalize 
the definition of the electron localization function (ELF) used to analyze bonding in 
ground- state systems (Becke 1990), to time-dependent processes: 

1 

TDELF (r, t) = (4.111) 

l + [D a (r,t)/D°(r,t)] Z 

with 

r\ / f \ ( f \ l |VMr,Q| 2 jl(rj) 

D a (r,t) = T a (r,t) - (4.112) 

4 n a (r,t) n a (r,t) 

where j a is the magnitude of the KS current-density of spin a, and x a (r, t) = 
Xfii \V<Pia( r i t)\ 2 is the KS kinetic energy-density of spin a. In Eq. 4.111, 
D®(r, t) = T% om (n a (r, t)) = |(6tt 2 ) 2/3 ^ /3 is the kinetic energy-density of the 
uniform electron gas. Using ALDA to evaluate the TDELF, this function has been 
useful for understanding time-resolved dynamics of chemical bonds in scattering 
and excitation processes (Burnus et al. 2005; Castro et al. 2007); features such as 
the temporal order of processes, and their time scales are revealed. As an example, 
in Fig. 4.5 we reproduce snapshots of the TDELF for laser-induced excitation of the 
7T — > 7r* transition in the acetylene molecule, studied in Burnus et al. (2005). Many 
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Fig. 4.5 Snapshots of the TDELF for the excitation of acetylene by a 17.5 eV laser pulse (from 
Burnus (2005)), polarized along the molecular axis. The pulse had a total length of 7 fs, an intensity 
of 1.2 xl0 14 W/cm 2 



interesting features can be observed from the TDELF. As the intensity of the laser 
field increases, the system begins to oscillate, and ionization is visible in the time 
slices at 1 .21 1 1 and 1 .5692 fs. The figure clearly shows that after 3 .5 fs, the transition 
from the ground- state to the antibonding state is complete: the original single torus 
signifying the triple bond in the ground-state has split into two separate tori, each 
around one carbon atom. 

General success for dynamics in strong fields has been slower than for linear 
response applications. There are three main reasons. First, many of the observables 
of interest are not simply related to the time-dependent one-body density, so that, 
in addition to the approximation for the xc functional, a new ingredient is needed: 
approximate "observable functionals" to extract the properties of interest from the 
KS system. Sometimes these are simply the usual quantum mechanical operators 
acting directly on the KS system, e.g. high-harmonic generation spectra are measured 
by the dipole moment of the system, J d 3 rn(r, co)z. But if the observable is not 
simply related to the density, such as ionization probabilities (Ullrich and Gross 
1997; Petersilka and Gross 1999), or cross-sections in atomic collisions (Henkel et 
al. 2009), in principle an observable functional is needed. Simply extracting double- 
ionization probabilities and momentum-densities using the usual operators acting on 
the KS wavefunction typically fails (Lappas and Van Leeuwen 1998; Wilken and 
Bauer 2006; Wilken and Bauer 2007; Rajam et al. 2009). 

Second, lack of memory dependence in the usual xc approximations has been 
suggested to be often far more problematic than in the linear-response regime as is 
discussed in Chap. 8. We must deal with the full xc potential fxcD?; @o]( r , 0» 
instead of the simpler xc kernel. The exact functional depends on the history of 
the density as well as on the initial state but almost all functionals used today 
are adiabatic. Third, a particularly severe difficulty is encountered when a system 
starting in a wavefunction dominated by a single Slater determinant evolves to a state 
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that fundamentally needs at least two Slater determinants to describe it. This is the 
time-dependent (TD) analog of ground-state static correlation, and arises in electronic 
quantum control problems (Maitra et al. 2002b; Burke et al. 2005a), in ionization 
(Rajam 2009), and in coupled electron-ion dynamics (Levine et al. 2006). The TD 
KS system evolves the occupied orbitals under a one-body Hamiltonian, remaining 
in a single Slater determinant: the KS one-body density matrix is always idempotent 
(even with exact functionals), while, in contrast, that of the true system develops 
eigenvalues (natural occupation numbers) far from zero or one in these applications 
(Appel and Gross 2010). The exact xc potential and observable functionals conse- 
quently develop complicated structure that is difficult to capture in approximations. 
For example, in Rajam et al.(2009), a simple model of ionization in two-electron 
systems showed that the momentum distribution computed directly from the exact 
KS system contains spurious oscillations due to using a single, necessarily delo- 
calized orbital, a non-classical description of the essentially classical two-electron 
dynamics. 



